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Asteroid  Orbit  Determination  and  Rotational  Period  Calculation  with  CCD  Astronomy 
Thesis  directed  by  Professor  Charles  Fosha,  Ph.D. 

This  paper  presents  data  collected  and  analyzed  relating  to  photometry  and 
astrometry  of  asteroids.  All  observations  were  accomplished  at  the  U.S.  Air  Force 
Academy  Observatory.  The  photometry  involves  determining  the  rotational  period  of 
asteroid  583  Klotilde.  Astrometry  was  performed  on  asteroid  1035  Amata  and  the 
calculated  position  was  used  to  determine  its  orbital  elements. 

Klotilde  was  selected  for  rotational  period  determination  based  on  its  relatively 
low  magnitude,  favorable  viewing  position,  and  no  previous  rotational  period 
information.  Two  hundred  six  images  of  Klotilde  were  taken  and  analyzed  over  four 
viewing  nights.  A  Photometries  (PM512)  Charge  Couple  Device  (CCD)  camera  attached 
to  a  61-cm  Cassegrain  telescope  was  used  for  these  observations.  Using  NOAO’s  IRAF 
software,  the  magnitudes  of  Klotilde  and  several  comparison  stars  were  determined. 
Using  an  Excel  spreadsheet,  differential  photometry  was  performed  and  the  lightcurve 
was  plotted.  The  four  nights  of  data  gave  a  9.210  ±  0.005  hour  synodic  period  with  an 
amplitude  of  0.18  magnitudes. 

Thirty-two  images  of  Amata  were  taken  on  six  different  viewing  nights.  The 
images  were  taken  with  an  ST-8  CCD  attached  to  a  41 -cm  Cassegrain  telescope.  The 
data  was  reduced  with  the  Astrometrica  software  package,  which  calculated  the  right 
ascension  (RA),  declination  (Dec),  and  magnitude  of  Amata  using  several  comparison 
stars.  The  computed  RA  and  Dec,  along  with  the  times  of  observation  were  then  used  to 


determine  the  orbital  elements  of  the  asteroid. 
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Gauss’  method  of  angles  only  orbit  determination  was  used  to  calculate  the  orbital 
elements.  Three  separate  computer  programs  were  used  to  calculate  the  orbital  elements. 
The  programs  used  were  the  Gauss-Encke-Merton  (GEM)  program  by  J.M.A.  Danby, 
ORBDET  by  Montenbruck  and  Pfleger,  and  FINDORB  by  Bill  Gray.  The  closest 
match  to  the  published  orbital  elements  came  form  FIND  ORB  which  used  a  batch  least 
squares  approach  to  use  all  of  the  observational  data.  ORBDET  and  GEM  used  only 
three  observations  at  a  time  and  yielded  less  accurate  results.  Due  to  suspected  errors  in 
the  GEM  program,  it  was  the  least  effective  method  used.  The  validity  of  our 
observations  was  verified  by  the  Minor  Planet  Center,  as  well  as  FIND  ORB  and 


ORBDET. 
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CHAPTER  I 


INTRODUCTION 


Astronomy  Background 

Our  solar  system  contains  thousands  of  objects  commonly  referred  to  as  asteroids 
or  minor  planets.  For  the  most  part,  they  are  just  rocks  in  space  that  orbit  around  the  Sun. 
There  is  much  theory  on  the  origin  of  asteroids,  from  the  break-up  of  planets,  to  pre- 
stellar  matter  that  failed  to  form  into  a  separate  planet.  Asteroids  can  be  found 
throughout  the  solar  system  in  a  variety  of  orbits.  The  vast  majority  are  located  in  fairly 
circular  orbits  located  between  the  orbits  of  Mars  and  Jupiter,  an  area  often  referred  to  as 
the  main  asteroid  belt.  Other  asteroids  of  particular  concern  have  orbits  that  intersect  that 
of  the  Earth,  which  are  referred  to  as  near  earth  asteroids.  Almost  all  asteroids  are 
relatively  small  with  diameters  of  less  than  five  miles  although  four  or  five  dozen  may 
exceed  one  hundred  miles. 1 

Throughout  this  thesis,  references  will  be  made  to  an  object’s  magnitude. 
Magnitudes  are  an  astronomer’s  method  of  indicating  brightness.  The  magnitude  scale  is 
an  inverse  scale  with  lower  numbers  indicating  brighter  objects.  Before  the  advent  of  the 
telescope,  all  observations  of  the  night  sky  were  obviously  made  with  only  the  naked  eye. 
Early  astronomers  classified  all  visible  stars  on  a  scale  of  1  to  6  based  on  their  apparent 
brightness.  Since  all  objects  are  different  distances  from  the  Earth  and  either  reflect  or 
generate  different  amounts  of  energy;  the  apparent  magnitude  is  not  an  accurate 
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estimation  of  its  true  luminosity.  For  some  observational  work  astronomers  correct 
apparent  magnitudes  to  absolute  magnitudes.  The  absolute  magnitude  is  the  apparent 
magnitude  a  star  would  have  if  it  were  at  a  distance  of  10  parsecs  or  32.6  light  years  from 
the  Earth.  All  of  the  magnitudes  discussed  in  this  thesis  will  refer  to  the  apparent 
magnitudes  measured  by  our  instruments. 

This  thesis  contains  two  distinct  sections.  The  first  section  deals  with  the 
photometry  of  an  asteroid  to  find  its  magnitude.  By  plotting  this  magnitude  as  a  function 
of  time,  a  rotational  period  for  the  asteroid  can  be  determined.  The  second  section 
describes  the  astrometry  of  a  different  asteroid  to  determine  its  right  ascension  and 
declination.  Using  these  coordinate  positions  and  the  observation  times,  Gauss’  method 
of  angles  only  orbit  determination  was  used  to  calculate  the  orbital  elements  of  the 
asteroid.  The  specifics  of  the  data  collection  procedures  and  data  analysis  will  be 
included  in  their  respective  chapters. 

Observational  Equipment 

All  observations  were  taken  at  the  U.S.  Air  Force  Academy  Observatory,  which 
can  be  seen  in  Fig.  1.1.  The  location  of  the  observatory  is  104°  52’  51.9”  W,  39°  00’ 
25.3”  N,  at  an  altitude  of  2180m.  A  large  portion  of  observational  astronomy  is 
conducted  with  Charge  Couple  Devices  (CCDs).  These  cameras  collect  photons  on  an 
electronic  chip  during  an  exposure.  After  the  shutter  closes  an  electric  potential  is  used  to 
“roll”  the  exposed  pixels  off  the  chip  so  that  they  can  be  counted.  CCD  astronomy  is 
considered  visual  astronomy  like  naked  eye  observations  because  it  makes  observations 
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in  the  visible  portion  of  the  spectrum.  Infrared  and  radio  are  two  types  of  astronomy  that 
observe  in  other  than  the  visible  spectrum. 


Fig.  1.1  USAF  Academy  Observatory 


CCD  photometry  was  accomplished  with  a  Photometries  (PM512)  camera 
attached  to  a  61 -cm  Cassegrain  telescope  with  an  equatorial  mount.  CCD  astrometry  was 
accomplished  with  a  41-cm  Cassegrain  telescope  with  an  equatorial  fork  mount.  The 
Cassegrain  configuration  indicates  that  the  telescope  achieves  a  longer  focal  length  by 
means  of  two  interior  mirrors.  An  illustration  of  a  Cassegrain  configuration  can  be  found 
in  Fig.  1.2.  An  equatorial  mount  aligns  the  prime  axis  of  the  telescope  with  the  Celestial 
North  Pole.  Therefore,  to  track  an  object  the  telescope  only  needs  to  rotate  in  one  axis. 
A  fork  mount  resembles  a  tuning  fork.  Depending  on  the  instrument  sizes,  there  can  be 
problems  viewing  objects  near  the  Celestial  North  Pole  with  a  fork  mount.  Pictures  of 
the  61-cm  and  41-cm  telescopes  can  be  seen  in  Fig.  1.3  and  1.4  respectively. 
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Fig.  1.2  Cassegrain  Configuration 
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Fig.  1.3  USAFA  61-cm  Telescope  Fig.  1.4  USAFA  41-cm  Telescope 


Photometry  vs.  Astrometry 

The  basic  principle  of  photometry  and  astrometry  are  the  same,  but  their  end  goals 
differ.  Both  processes  use  a  CCD  to  collect  visible  light  photons.  For  asteroids,  the  light 
is  generated  by  the  Sun  and  reflected  off  the  asteroid.  Neighboring  stars  obviously 
generate  their  own  light.  The  neighboring  stars  will  hereafter  be  referred  to  as 
comparison  or  reference  stars.  After  an  exposure  is  taken  the  photons  are  electronically 
counted  and  a  computer-generated  representation  of  the  objects  in  the  field  of  view  of  the 


camera  is  displayed. 
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In  photometry,  only  the  number  of  photons  and  therefore  the  magnitudes  of  the 
asteroid  and  comparison  stars  are  calculated.  The  data  are  then  differentially  corrected 
and  plotted.  The  astrometry  software  used  also  calculates  the  magnitude,  but  this  is  only 
a  secondary  calculation  as  far  as  our  research  was  concerned.  Our  main  purpose  was  to 
calculate  the  position  of  the  asteroid.  Since  all  stars  can  be  considered  fixed  in  space,  the 
astrometry  software  can  use  the  known  coordinates  of  several  of  the  comparison  stars  to 
determine  the  coordinates  of  the  asteroid  in  question. 

The  calculated  positions  are  given  in  units  of  right  ascension  and  declination. 
Right  ascension  and  declination  are  the  same  for  any  object  regardless  of  viewing  time  or 
viewing  location.  This  feature  makes  them  the  units  of  choice  for  astronomers  who  can 
use  the  coordinates  to  determine  when  a  particular  object  is  visible  from  their  observing 


location  if  it  is  at  all. 


CHAPTER  II 


583  KLOTILDE  PHOTOMETRY  OBSERVATION  PROCEDURES 

Selection  of  583  Klotilde 

The  purpose  of  the  photometry  research  was  to  select  an  asteroid  to  observe 
whose  rotational  period  had  not  previously  been  determined.  To  begin,  Guide  software 
version  3.0  was  used  to  identify  asteroids  that  would  be  visible  from  our  location.  The 
Guide  software  incorporates  data  from  the  Hubble  Guide  Star  Catalog.2  Although  these 
programs  have  different  names,  they  are  quite  different.  In  this  thesis.  Guide  will  refer  to 
the  program  used  to  see  the  positions  of  asteroids  and  stars  at  different  viewing  times. 
The  Hubble  Guide  Star  Catalog  is  a  database  of  celestial  objects  used  by  Guide  as  well  as 
other  programs  mentioned  later.  Guide  allowed  us  to  input  the  viewing  time  and  location 
so  we  could  see  which  asteroids  were  visible  during  our  observation  time.  Obviously  the 
asteroids  we  were  considering  had  already  been  discovered  but  due  to  the  large  number 
of  asteroids  in  our  solar  system,  several  of  them  had  not  been  observed  extensively.  The 
Guide  Star  Catalog  also  contains  the  magnitudes  for  the  asteroids  we  were  considering. 
Since  a  higher  magnitude  star  is  dimmer,  we  arbitrarily  set  an  upper  magnitude  limit  of 
about  19,  which  would  be  close  to  the  limit  at  which  our  telescope  could  detect  the 
asteroid. 

From  a  list  of  approximately  five  candidate  stars  we  consulted  information  by 
Harris  to  see  which  of  those  asteroids  did  not  already  have  a  rotational  period 


7 


determined.3  Based  on  this  information,  asteroid  583  Klotilde  was  selected.  Klotilde  was 
discovered  on  1898  Dec  31  in  Vienna  by  J.  Palisa.4  Its  approximate  diameter  can  be 
estimated  at  about  86  km  based  on  its  magnitude;  however,  no  guess  as  to  its  shape  can 
be  made.5 

Planning  the  Observations 

Once  538  Klotilde  was  selected,  we  needed  to  plan  the  observations.  Again 
consulting  Guide,  we  were  able  to  adjust  the  time  settings  for  the  night  we  planned  to 
observe.  The  software  showed  us  the  asteroid  with  respect  to  the  other  stars  that  would 
be  in  the  field  of  view  during  that  time.  Since  the  CCD  has  a  field  of  view  of  only  a  few 
arcminutes,  it  was  important  to  ensure  that  there  was  at  least  one  suitable  comparison  star 
in  the  field  of  view  with  the  asteroid.  A  suitable  comparison  star  should  be  relatively 
bright  with  a  magnitude  of  approximately  10  to  17.  One  comparison  star  is  essential; 
however,  two  or  three  are  desirable  as  will  be  seen  later  in  the  discussion  of  data 
reduction.  If  there  weren’t  any  suitable  comparison  stars  on  a  particular  night  then 
changing  observing  nights  was  considered. 

Another  factor  to  consider  is  that  over  the  course  of  the  night  the  asteroid  will  be 
moving.  By  changing  the  settings  in  Guide  to  several  hours  later,  we  could  see  the 
general  direction  and  distance  it  would  travel.  It  was  necessary  to  make  sure  that  the 
asteroid  didn’t  travel  too  far  away  from  the  comparison  stars  so  that  they  would  no  longer 
be  in  the  same  field  of  view.  We  also  checked  to  make  sure  the  asteroid  would  not  pass 
too  close  to  another  star  which  would  interfere  with  determining  which  photons  came 
from  which  object.  All  of  these  considerations  needed  to  be  made  along  with  the  weather 
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forecast.  Observing  on  consecutive  nights  would  make  determining  the  rotational  period 
easier,  particularly  for  asteroids  with  longer  periods.  Therefore  selecting  an  observing 
night  was  a  balance  between  favorable  comparison  stars  on  consecutive  nights  and  the 
weather.  Invariably  something  unexpected  will  occur  while  observing;  therefore,  backup 
nights  are  always  a  good  idea. 

Observing  Klotilde 

A  61 -cm  Cassegrain  telescope  was  used  for  all  photometry  observations.  Prior  to 
observations  each  night,  the  CCD  needed  to  be  cooled  with  liquid  nitrogen.  A 
temperature  of  -120°  C  was  desired  for  observations.  The  telescope  was  also  focused  to 
give  a  clear  image.  All  of  the  images  taken  of  583  Klotilde  had  a  five-minute  exposure 
time  using  a  standard  Johnson  R  band  (red)  filter.  An  exposure  of  this  length  makes  any 
error  due  to  shutter  inefficiency  negligible.  The  images  were  taken  using  PMIS  software. 
Four  nights  of  observations  were  taken  of  Klotilde.  See  table  2.1  for  a  list  of  the 
observing  nights,  the  time  of  observations,  number  of  observations,  and  coordinates  of 
the  asteroid  that  night.  Note  that  the  right  ascension  and  declination  refer  to  the  telescope 
pointing  coordinates  not  the  coordinates  for  Klotilde  specifically.  Since  the  asteroid  is 
orbiting  the  Sun  and  not  fixed  in  space  like  stars  are,  these  pointing  coordinates  are 
actually  changing  continuously  over  the  four  or  five  hours  that  observations  were  made. 
Over  these  intervals  however,  the  changes  are  not  particularly  relevant. 
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UT  98  Feb  13 

UT  4.8 

UT  10.3 

48 

8  53  76.9 

5  49  46 

UT  98  Feb  22 

UT  3.4 

UT  9.0 

52 

8  46  36.7 

6  16  30 

UT  98  Mar  12 

UT  1.9 

UT  9.4 

56 

8  38  34.0 

7  07  16 

UT  98  Mar  13 

UT  2.8 

UT  7.5 

50 

8  38  06.1 

7  12  33 

Table  2.1.  Photometry  Observing  Log 


With  each  observation,  an  observing  log  was  filled  out  with  other  procedural 
information.  The  logs  from  all  observations  at  the  Air  Force  Academy  Observatory  are 
maintained  in  a  binder  for  possible  future  reference.  Included  in  the  log  is  the  start  time, 
read  to  the  nearest  second  (Universal  Time)  from  a  wristwatch  synchronized  with  the 
U.S.  Naval  Observatory  master  clock.  The  exposure  duration  is  specified  with  the 
camera  software.  The  pointing  coordinates  of  the  telescope  during  the  exposure  are  also 
recorded.  Although  the  telescope  control  software  automatically  tracks  the  same  position' 
in  the  sky  correcting  for  the  Earth’s  rotation,  small  offset  adjustments  are  needed 
periodically.  In  the  event  that  the  telescope  needs  to  be  moved  for  something  such  as 
viewing  another  object  or  refilling  the  liquid  nitrogen  dewar  it  could  be  repositioned  later. 
The  filter  used  was  also  recorded.  For  some  astronomy  purposes  different  filters  are  used 
for  the  same  image,  but  as  previously  mentioned,  all  of  our  exposures  used  the  R  filter. 

Other  useful  information  recorded  in  the  log  is  the  hour  angle,  which  indicates 
distance  west  (positive)  or  east  (negative)  from  the  meridian  passing  through  the  Celestial 
North  Pole  and  the  local  zenith  position.  It  is  not  used  directly  in  data  analysis,  but  can 
be  useful  in  determining  how  late  observations  can  be  made.  Also  recorded  was  the 
airmass,  which  indicates  the  amount  of  atmosphere  through  which  the  telescope  is 
looking.  Looking  directly  overhead,  the  telescope  is  looking  through  one  airmass.  For 
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angles  off  zenith,  the  amount  of  atmosphere  the  telescope  is  looking  through  increases. 
The  airmass  also  increases  accordingly.  As  with  the  hour  angle,  the  airmass  is  not  used 
directly  in  data  analysis  but  can  be  useful  as  a  reference.  If  an  image  is  blurred  or  hazy 
the  airmass  can  be  checked  to  see  how  close  the  object  was  to  the  horizon  at  the  time  of 
the  exposure.  Light  from  objects  closer  to  the  horizon  has  more  atmosphere  to  pass 
through  and  can  be  more  distorted. 

In  addition  to  the  actual  images  of  Klotilde  approximately  five  flat  field  images 
and  ten  bias  frames  were  taken  each  night.  Flat  field  corrections  take  into  account  the 
pixel-to-pixel  sensitivity  as  well  as  other  unavoidable  effects  such  as  vignetting  or  dust 
on  the  CCD.6  The  flat  field  images  were  25-second  exposures  taken  either  of  the 
illuminated  inside  of  the  dome  or  the  twilight  sky.  Since  new  flat  field  images  are  needed 
each  night  observations  are  made,  both  of  these  techniques  were  used  depending  on  the 
time  of  observations.  The  bias  frames  are  zero  second  dark  current  exposures.  A  similar 
exposure  called  a  dark  frame  (discussed  in  the  astrometry  section,  Chapter  V)  was  not 
needed  since  the  CCD  was  cooled  to  less  than  -100°  C.  The  purpose  of  these  additional 
frames  will  be  discussed  along  with  data  reduction. 


CHAPTER  III 


583  KLOTILDE  DATA  REDUCTION 

Data  Reduction  Software 

The  PMIS  software  saves  the  CCD  images  in  the  Flexible  Image  Transport 
System  (FITS)  format7,  which  is  the  accepted  format  for  professional  astronomers.6 
These  files  are  then  transferred  to  a  Sun  computer  workstation  where  the  data  reduction  is 
actually  accomplished.  The  software  used  to  reduce  the  images  is  DIGIPHOT  8’ 9,  which 
is  also  the  standard  among  professional  astronomers.6  Specifically  we  used  aperture 
photometry,  which  is  a  package  in  the  DIGIPHOT  software.  This  software  runs  on  the 
Image  Reduction  and  Analysis  Facility  (IRAF)  package  available  through  the  National 
Optical  Astronomy  Observatories  (NOAO). 

Image  Corrections  and  Analysis 

After  importing  the  FITS  images  and  organizing  them  into  user  specified 
directories,  they  must  be  converted  to  the  proper  image  format  used  by  IRAF.  The  bias 
and  flat  field  images  as  well  as  the  light  images  are  all  converted  to  this  format.  Next,  a 
master  bias  image  is  created  from  an  average  of  all  the  bias  frames.  IRAF  has  a 
command  that  allows  the  user  to  define  the  parameters  and  then  it  automatically  performs 
this  operation.  This  master  bias  is  then  subtracted  from  each  of  the  light  images  as  well 
as  the  flat  field  images.  A  function  similar  to  the  one  used  to  create  the  master  bias 
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image  is  used  to  create  a  master  flat  field  image.  It  is  then  necessary  to  normalize  the 
average  pixel  value  in  the  master  flat  field  to  1  before  applying  the  flat  field  correction 
into  each  of  the  light  images.  The  light  images  are  now  ready  to  be  displayed  for 
analysis. 

A  picture  of  the  image  is  displayed  on  the  screen.  After  viewing  a  few  images 
from  the  beginning  and  end  of  the  night  it  is  possible  to  determine  which  stars  will  be 
used  as  comparison  stars.  The  comparison  stars  should  be  similar  in  magnitude  (visual 
inspection  is  sufficient),  and  should  appear  in  all  or  most  of  the  frames.  It  is  also 
important  to  select  comparison  stars  that  are  not  too  close  to  the  edges  of  the  frame  since 
this  can  affect  the  pixel  count.  For  our  purposes,  one  to  three  comparison  stars  were 
selected  for  each  observing  night.  An  example  of  a  raw  image  from  the  PMIS  software 
can  be  seen  in  Fig.  3.1.  Notice  the  donut  shaped  smudges  in  parts  of  the  image.  These 
smudges  are  removed  by  applying  the  flat  field  correction.  Fig.  3.2  contains  the  same 
image,  which  has  been  enhanced  to  simulate  the  bias  and  flat  field  corrections  by 
adjusting  the  contrast  scale.  Klotilde  is  the  object  indicated  by  the  yellow  arrow.  The 
two  other  stars  identified  with  red  arrows  were  used  as  comparison  stars.  The  bright  star 
towards  the  top  indicated  by  the  blue  arrow  could  also  have  been  used,  but  the  reason  it 


was  not  will  be  discussed  later. 


Fig.  3.2  Enhanced  PMIS  583  Klotilde  Image.  FOV:  3.7  X  3.7  Arcminutes 


14 


The  images  then  need  to  be  inspected  for  cosmic  rays.  One  of  the  disadvantages 
of  CCD  photometry  with  long  exposures  is  that  it  is  susceptible  to  pollution  from  cosmic 
or  galactic  rays.  The  cosmic  rays  usually  appear  as  highly  concentrated  bright  pixels, 
which  stand  out  from  the  surrounding  dark  background.  Some  cosmic  ray  spikes  can  be 
seen  as  single  bright  pixels  in  Fig  3.2.  The  stars  and  asteroids  are  quite  large  compared  to 
the  cosmic  ray  spikes  making  them  easy  to  identify.  Only  those  spikes  in  close  proximity 
to  the  asteroid  and  comparison  stars  are  of  concern.  Close  proximity  is  difficult  to  define 
but  a  sense  of  this  is  developed  after  processing  a  few  images.  If  the  spikes  are  located 
extremely  close  to  the  star  or  asteroid  however,  it  can  be  difficult  to  distinguish  the  spike 
from  the  object’s  photons.  Since  the  cosmic  rays  will  increase  the  photon  count  for  the 
object  they  are  close  to,  it  is  necessary  to  eliminate  them  from  the  image.  The 
DIGIPHOT  software  allows  the  user  to  replace  the  pixels  affected  by  the  cosmic  rays 
with  the  average  value  of  the  surrounding  pixels.  This  is  only  an  approximation,  but  it 
yields  a  pixel  value  much  more  realistic  than  the  cosmic  ray  value. 

After  the  cosmic  ray  spikes  of  concern  have  been  eliminated,  the  user  identifies 
the  center  of  the  object  desired  with  the  mouse  and  an  aperture  around  the  object  is 
specified  to  count  the  photons  in.  DIGIPHOT  then  counts  the  total  number  of  photons 
within  the  specified  aperture.  Usually  several  aperture  sizes  are  selected  and  the  best  one 
is  selected  later.  A  plot  of  the  photons  is  then  displayed.  This  is  inspected  to  ensure  no 
remaining  cosmic  ray  spikes  within  the  specified  apertures  are  corrupting  the  photon 
count.  The  cosmic  rays  will  lie  outside  the  main  photon  curve  since  they  have  a  much 
higher  intensity  than  the  surrounding  pixels.  If  additional  cosmic  ray  spikes  are  detected 
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they  need  to  be  eliminated  as  discussed  previously  and  the  photometry  process  is 
repeated. 

In  each  frame,  photometry  needs  to  be  accomplished  for  the  asteroid  and  each 
comparison  star.  It  is  important  to  select  the  asteroid  and  comparison  stars  in  the  same 
order  for  each  frame.  The  result  of  this  process  is  a  file,  which  includes  the  calculated 
magnitude  of  each  object  selected,  and  the  error  for  each  calculation.  The  magnitude  and 
error  based  on  each  aperture  specified  is  calculated.  Usually  the  calculations  from  a 
larger  aperture  are  used  unless  the  selected  object  is  too  close  to  the  edge  of  the  frame.  In 
this  case,  a  smaller  aperture  size  may  be  used.  A  little  experimentation  provides  the 
proper  balance  of  including  all  the  photons  in  the  aperture  while  still  remaining  in  the 
frame  boundaries. 

From  this  file,  the  magnitude  and  error  for  each  object  are  recorded  in  a  notebook. 
This  data  is  then  used  to  perform  the  differential  photometry  analysis  to  determine  the 
rotational  period. 


CHAPTER  IV 


583  KLOTILDE  DATA  ROTATIONAL  PERIOD  DETERMINATION 

Lightcurve  Information 

The  lightcurve  referred  to  in  this  thesis  is  a  plot  of  the  magnitude  versus  time. 
The  purpose  for  plotting  the  lightcurve  of  an  asteroid  is  based  on  the  fact  that  asteroids 
have  irregular  shapes.  For  the  purpose  of  illustration  assume  that  an  asteroid  has  a 
roughly  potato  shape.  As  the  asteroid  rotates,  different  faces  or  sides  of  it  face  the  Earth. 
When  an  edge  with  a  larger  cross  sectional  area  faces  the  Earth,  it  will  reflect  more  light 
than  an  end  with  a  smaller  cross  sectional  area.  When  more  photons  are  reflected 
towards  the  Earth  the  CCD  registers  a  higher  photon  count  and  thus  a  lower  (brighter) 
magnitude  is  determined  for  that  asteroid.  As  the  asteroid  makes  one  complete  rotation, 
the  sides  facing  the  Earth  will  be  edge,  end,  edge,  end. 

The  potato  asteroid  example  assumes  that  the  two  edges,  or  put  another  way,  the 
front  and  the  back,  have  roughly  the  same  cross  sectional  area,  as  do  the  two  ends.  Under 
these  conditions,  a  plot  of  the  magnitudes  should  reveal  a  figure  similar  to  a  sine  wave 
with  two  maxima  and  two  minima.  By  matching  up  the  location  where  the  curve  begins 
to  repeat  itself,  the  rotational  period  can  be  determined.  If  the  rotational  period  is  such 
that  most  of  a  complete  rotation  can  be  observed  in  one  or  two  viewing  nights  it  should 
be  possible  to  determine  the  rotational  period.  In  Fig.  4.1,  the  period  (where  the  curve 
repeats  itself)  is  the  entire  length  of  the  frame, 
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Hagnitude 


Time 

Fig.  4.1  Potato  Shaped  Asteroid  Lightcurve 

Of  course  in  reality,  most  asteroids  are  not  as  symmetric  as  the  potato  example.  If 
one  edge  or  end  has  a  slightly  larger  cross  sectional  area  than  the  other  does,  then  the 
light  reflected  and  therefore  the  respective  magnitudes  will  be  different.  In  this  case,  one 
of  the  maximums  or  minimums  might  be  larger  than  the  other,  as  shown  in  Fig.  4.2. 
However,  in  most  circumstances  it  is  still  reasonable  to  assume  that  a  double  maxima  and 
minima  lightcurve  will  result. 


Magnitude 


Time 

Fig.  4.2  Irregularly  Shaped  Asteroid  Lightcurve 

Under  certain  circumstances,  it  is  impossible  to  determine  a  rotational  period  with 
any  certainty.  If  the  asteroid  should  happen  to  be  roughly  spherical  like  a  planet, 
(improbable  since  many  asteroids  are  loose  conglomerations  of  smaller  chunks  and  rocks, 
or  the  result  of  catastrophic  collisions)  then  relatively  no  change  would  be  detected  in  the 
plotted  magnitude  making  a  rotational  period  determination  extremely  difficult.  Also,  if 
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the  asteroid  is  rotating  extremely  slowly,  then  several  nights  of  observations  are 
necessary  just  to  observe  a  fraction  of  the  period.  This  would  make  it  a  poor  candidate 
for  observation.  The  main  reason  for  this  is  that  several  consecutive  nights  of 
observations  would  be  necessary  which  is  complicated  by  trying  to  find  consecutively 
clear  nights.  After  plotting  one  night’s  data,  it  would  be  possible  to  conclude  whether  or 
not  this  is  a  good  asteroid  for  further  observations.  The  results  from  such  a  scenario  could 
still  be  useful  however,  by  simply  stating  that  the  named  asteroid  has  a  rotational  period 
of  at  least  X  hours.  This  could  be  of  benefit  to  a  future  astronomer  who  is  trying  to  find 
an  asteroid  with  a  relatively  short  period. 

Differential  Photometry  Purpose  and  Calculations 

The  magnitude  and  error  data  recorded  from  the  DIGIPHOT  calculations  were 
then  fed  into  an  Excel  worksheet.  The  worksheet  for  each  observing  night  can  be  found 
in  Appendix  A.  The  format  for  this  worksheet  was  written  by  Maj.  Charles  Wetterer 
from  the  Air  Force  Academy  Department  of  Physics  for  use  in  differential  photometry. 
Differential  photometry  is  the  process  of  determining  a  difference  in  magnitudes  based  on 
observational  calculations.  The  reason  this  process  is  needed  is  that  the  observations  are 
made  over  a  long  period  of  time  (four  or  five  hours  for  our  data)  for  each  observing  night. 
Over  this  period,  it  is  likely  that  clouds  will  drift  in  and  out  of  the  field  of  view  of  the 
telescope,  which  will  reduce  the  number  of  photons  reaching  the  CCD.  There  is  also  a 
change  in  the  airmass  over  the  viewing  period,  which  affects  the  number  of  photons 
collected.  Differential  photometry  considers  these  factors  and  adjusts  the  calculated 
magnitude  accordingly. 
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The  differential  corrections  are  applied  to  each  object  (asteroid  and  comparison 
stars)  against  each  other  object.  First,  the  magnitude  of  each  comparison  star  is 
subtracted  from  the  asteroid 


mK-  nn 

(4.1) 

mK  -  m2  -(mi  -  m2) 

(4.2) 

mK  -  m3  -  (nii-m3) 

(4.3) 

where  K  represents  Klotilde  and  the  numbers  represent  the 
comparison  stars 

The  third  term  subtracted  from  equations  4.2  and  4.3  is  used  to  adjust  the  values  to  the 
first  comparison  star.  This  term  is  actually  the  average  of  all  of  the  individual 
observations.  The  time  is  computed  at  the  middle  of  the  exposure  time  as 

Hour  +  (Min  /  60)  +  (Sec  /  3600)  +  (Exposure  length  /  (3600  *  2))  (4.4) 

An  offset  is  applied  to  correct  for  the  change  in  observed  magnitudes  over  the  entire 
night.  The  offset  is  the  magnitude  of  the  first  comparison  star,  which  is  taken  from  the 
Guide  Star  Catalog.  The  offset  is  added  to  equation  4.1  to  get  the  combined  magnitude. 


mComb  =  (rnK  -  mO  +  offset 


(4.5) 
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The  combined  error  is  the  square  root  of  the  sum  of  the  squares  of  the  individual  errors, 


which  were  calculated  by  DIGIPHOT. 


CTcomb  =  Sqrt  (a2K  +  ct2i) 


(4.6) 


The  combined  magnitude  of  the  asteroid  with  the  offset  correction  is  then  plotted  versus 


time  to  obtain  the  lightcurve. 


Check  for  Variable  Stars 

The  magnitude  of  the  comparison  stars  are  also  corrected  against  each  other  and 
plotted  to  determine  if  any  of  them  are  variable.  An  example  of  such  a  plot  is  shown  in 
Fig.  4.3.  The  plot  should  be  roughly  a  horizontal  line.  Should  there  be  fluctuations  in  the 
magnitude,  it  is  possible  that  one  of  the  stars  is  variable.  A  variable  star  is  one  whose 
light  output  varies  due  to  a  variety  of  reasons  that  will  not  be  discussed  here. 


Comparison  Star  Magnitude  Plot 
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Fig.  4.3  Comparison  Star  Magnitude  Plot 
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It  is  interesting  to  note  that  the  data  from  one  night’s  observations  indicated  that 
one  of  the  comparison  stars  was  indeed  a  variable  star.  This  is  the  star  in  the  upper  part 
of  Fig.  3.2  indicated  by  the  blue  arrow.  See  Fig.  4.4  for  the  comparison  star  magnitude 
plot  for  that  star.  The  star  was  checked  against  a  published  list  of  variable  stars,  but  was 
not  found  on  it.  With  the  suspicion  of  being  variable,  that  star  was  not  used  as  a 
comparison  star  for  differential  photometry.  Luckily,  there  were  two  other  candidates  for 
comparison  stars  from  that  night’s  data.  Determining  the  exact  nature  of  the  variability 
of  that  star  provides  and  excellent  opportunity  for  further  astronomy  research.  This  star’s 
designation  is  GSC  223:1761  and  is  located  in  the  constellation  Cancer. 


Variable  Comparison  Star  Magnitude  Plot 
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Fig.  4.4  Variable  Comparison  Star  Magnitude  Plot 
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Analyzing  Klotilde’s  Lightcurve 

After  the  differential  corrections  were  made,  the  magnitude  of  the  asteroid  was 
plotted  versus  time.  The  data  from  the  first  night  indicated  what  appeared  to  be  a  double 
minima  and  a  double  maxima.  We  were  concerned  however,  that  the  amplitude  was  only 
about  0.18  magnitudes.  In  addition,  the  comparison  star  plot  did  not  resemble  what  we 
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expected  it  to  look  like.  We  decided  to  attempt  another  night  of  observations  on  Klotilde 
to  see  if  new  data  could  confirm  that  the  data  from  the  first  night  was  reliable.  The 
second  night’s  data  indicated  a  very  distinct  and  smooth  curve  varying  only  towards  the 
end  of  the  evening  as  Klotilde  was  moving  close  to  the  horizon.  This  data  also  produced 
a  double  minima  and  double  maxima  lightcurve.  The  comparison  star  plot  was  also 
much  more  acceptable. 

With  this  information,  it  appeared  as  if  Klotilde  possessed  a  period  of 
approximately  4.5  hours.  This  was  difficult  to  state  with  certainty,  however,  since  more 
than  a  week  had  passed  between  the  two  nights  of  observations.  This  allowed  for  other 
integer  multiples  of  the  period  based  on  the  number  of  aliases.  The  number  of  aliases 
results  from  the  frequency  of  our  sampling  against  the  entire  period.  If  the  observation 
rate  is  too  small  then  the  period  indicated  by  our  data  may  not  cover  one  entire  period. 
Therefore,  two  more  consecutive  nights  of  observations  were  conducted.  Observing  on 
consecutive  nights  reduces  the  number  of  possible  aliases. 

The  third  night  of  observations  was  our  longest  observing  night.  The  data  from 
that  night  indicated  a  triple  minima  and  maxima  lightcurve.  The  period  associated  with 
this  lightcurve  was  approximately  9.2  hours.  A  triple  minima  and  maxima  lightcurve  is 
unusual,  which  made  us  suspicious  because  an  asteroid  with  such  a  lightcurve  would 
have  a  particularly  unusual  shape.  Adding  the  data  from  all  four  nights  confirmed  our 
data  from  the  third  night  that  the  lightcurve  did  contain  a  triple  minima  and  maxima.  By 
adjusting  the  period  offset,  we  were  able  to  find  the  period  by  visually  determining  when 
the  curves  from  each  night  best  lined  up.  The  most  accurate  period  we  could  determine 
from  our  data  was  9.210  +  0.005  hours.  This  is  a  synodic  not  sidereal  period.  A  sidereal 
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period  is  the  time  it  takes  for  the  asteroid  to  return  to  a  reference  position  in  inertial 
space.  The  synodic  period  is  the  time  for  the  same  part  to  the  asteroid  to  be  facing  the 
Earth.  Since  we  are  observing  from  the  Earth,  this  is  the  type  of  period  we  are  able  to 
calculate.  The  synodic  period  may  be  either  shorter  or  longer  than  the  sidereal  period 
depending  on  the  direction  of  rotation  of  the  asteroid  and  its  position  in  the  solar  system 
relative  to  the  Earth.  The  amplitude  without  error  considered  was  determined  to  be  0.18 
magnitudes.  A  graph  of  the  final  lightcurve  can  be  found  in  Fig.  4.5. 


Fig.  4.5  Final  Lightcurve  Plot  for  583  Klotilde 


Results  Submitted  for  Publication 

Since  the  rotational  period  of  583  Klotilde  had  not  previously  been  recorded,  we 
submitted  an  article  summarizing  our  findings  to  the  Association  of  Lunar  and  Planetary 
Observers  for  publication  in  their  quarterly  newsletter.  The  Minor  Planet  Bulletin. 
Should  our  article  be  accepted  for  publication,  it  would  most  likely  be  published  in 
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August.  Co-authoring  the  article  was  Maj.  Charles  Wetterer.  A  copy  of  our  submitted 
paper  can  be  found  in  Appendix  B. 


CHAPTER  V 


1035  AMATA  ASTROMETRY  OBSERVATION  PROCEDURES 

Selecting  1035  Amata 

CCD  Astrometry  was  conducted  to  determine  the  right  ascension  and  declination 
of  several  asteroids.  The  first  part  of  this  process  was  to  select  an  asteroid  or  asteroids  on 
which  to  conduct  astrometry.  It  was  decided  to  use  the  critical  list  of  asteroids  published 
by  the  Minor  Planet  Center  to  determine  which  asteroids  to  observe.  The  critical  list  is  a 
regularly  updated  publication,  which  is  used  by  astronomers  to  select  asteroids  for 
observing.  The  basic  reason  that  asteroids  are  put  on  the  critical  list  is  because  there  are 
relatively  few  observations  of  that  asteroid,  and  more  are  needed  to  get  an  accurate 
estimation  of  the  orbit.  Specifically,  there  are  six  different  categories  of  asteroids  on  the 
critical  list.10  The  categories  are  broken  down  based  on  the  number  of  oppositions  the 
asteroid  has  been  observed  at  over  a  certain  time  frame.  The  specific  asteroid  that  we 
eventually  selected  was  1035  Amata,  although  it  wasn’t  selected  until  later.  Amata  was 
listed  in  Category  5  of  the  critical  list  which  is  for  objects  observed  at  four  or  more 
oppositions  only  one  night  in  the  last  ten  years.10  Opposition  refers  to  the  asteroid’s 
location  in  its  orbit  with  respect  to  the  Earth  and  the  Sun.  Since  the  asteroid’s  orbit  lies 
outside  the  Earth’s,  opposition  occurs  when  the  Earth  and  asteroid  both  lie  on  the  same 
side  of  the  Sun.  Conjunction,  the  opposite  of  opposition,  occurs  when  the  Sun  lies 
between  the  Earth  and  the  asteroid.  See  Fig  5.1  for  an  illustration  of  this  concept. 
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There  were  two  reasons  that  we  decided  to  view  asteroids  from  the  critical  list. 
First,  it  could  be  assumed  that  the  critical  list  asteroids  would  be  in  favorable  viewing 
positions  during  our  observation  times.  In  addition  to  obtaining  data  useful  for  this 
thesis,  the  data  we  collected  could  be  reported  to  the  Minor  Planet  Center  for  use  in 
updating  their  asteroid  ephemerides.  Starting  with  the  critical  list  we  next  consulted 
Guide  to  see  when  specific  asteroids  would  be  visible  during  our  viewing  time.2  Guide 
allows  the  observer  to  select  the  asteroid  number,  which  it  then  shows  in  relation  to  other 
nearby  stars  at  the  time  specified.  It  was  possible  to  set  the  time  to  when  we  would  be 
observing  to  see  if  there  were  favorable  comparison  stars  in  close  vicinity  to  the  asteroid. 
At  least  four  comparison  stars  are  highly  recommended  and  seven  to  twelve  comparison 
stars  are  preferred  by  the  analysis  software  used  to  determine  the  position  of  the  asteroid. 

13 

Special  consideration  was  given  to  the  fact  that  the  asteroids  would  have  some 
relative  motion  over  the  period  that  we  would  be  observing  it.  It  was  therefore  necessary 
to  make  sure  that  the  asteroid  would  not  be  passing  too  close  to  other  stars  over  the 
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observing  time.  In  order  to  perform  the  astrometry  we  would  need  two  pairs  of 
observations  separated  by  approximately  one  hour.  The  reason  for  this  separation  is  so 
that  visual  inspection  could  ensure  that  the  object  suspected  to  be  the  asteroid  was  in  fact 
our  target  asteroid.  Since  the  asteroid  will  have  relative  motion  over  one  hour  and  stars 
do  not,  the  object  that  has  moved  from  one  observing  time  to  the  next  is  the  target 
asteroid. 

Another  area  of  consideration  was  the  magnitude  of  the  asteroid.  Most  asteroids, 
because  of  their  distance  from  the  Sun  and  size,  will  have  magnitudes  in  the  range  of  13 
to  22.  Due  to  the  limitations  of  our  telescope,  we  decided  a  reasonable  upper  magnitude 
limit  would  be  about  19  magnitudes.  The  problem  with  observing  dimmer  asteroids  is 
that  it  might  be  difficult  to  get  an  accurate  image  with  an  exposure  time  of  only  a  few 
minutes. 

Initially  several  asteroids  were  selected  for  observing.  This  was  done  to  provide 
backup  data  should  the  data  on  any  specific  asteroid  not  be  as  good  as  another  was.  This 
allowed  us  to  wait  until  later  to  select  the  specific  asteroid  that  we  would  use  for  orbit 
determination.  Additionally  any  observations  we  could  make  would  be  useful  to  the 
Minor  Planet  Center.  Since  there  needed  to  be  an  hour  separation  between  pairs  of 
observations  for  each  asteroid,  we  could  make  better  use  of  that  time  by  imaging  other 
asteroids. 

Observing  Asteroids 

The  observations  were  taken  using  a  41-cm  Cassegrain  telescope  and  an  ST-8 
CCD  camera.  This  CCD  was  electronically  cooled  to  -10°C.  The  temperature  was 
specified  with  Sky-Pro  CCD  Astronomy  Software  which  is  also  the  software  that  took  the 
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exposures.1 1  After  the  telescope  was  initialized  to  know  what  coordinates  it  was  pointing 
at,  it  was  necessary  to  focus  the  telescope.  The  Sky-Pro  software  has  a  focusing  option  to 
assist  in  obtaining  a  good  focus.  A  random  star  is  chosen  and  highlighted  on  the 
computer  screen.  Short,  repeated  exposures  are  taken  automatically  by  the  camera  while 
the  operator  adjusts  the  focus.  By  visual  inspection  of  the  images  resulting  from  the 
different  focuses,  the  operator  can  select  the  best  one. 

To  initially  locate  the  asteroid  the  telescope  was  slewed  to  the  right  ascension  and 
declination  given  by  Guide.  An  exposure  just  long  enough  to  pick  up  the  stars  in  the  field 
of  view  (approximately  20  seconds)  was  taken.  This  image  was  then  compared  against 
the  picture  given  by  Guide.  Some  offset  corrections  were  necessary  to  center  the  asteroid 
and  the  planned  comparison  stars  in  the  field  of  view. 

For  most  of  the  asteroids  observed  a  two-minute  exposure  time  was  sufficient. 
For  one  asteroid  with  a  higher  magnitude,  a  three-minute  exposure  was  taken.  For  each 
observation  taken  there  were  two  frames  that  were  exposed,  a  light  frame  and  a  dark 
frame.  A  light  frame  is  obtained  by  actually  opening  the  camera  shutter  and  collecting 
photons  from  the  sky.  The  dark  frame  is  an  exposure  of  the  same  length  as  the  light 
frame,  but  is  taken  with  the  shutter  closed.  The  dark  frame  is  necessary  since  the  camera 
is  not  operating  below  -100°C.6  It  is  needed  to  determine  the  amount  of  dark  current  or 
thermal  noise  in  the  system  and  the  DC  offset  voltage  (bias)  applied  to  the  CCD,  which  is 
later  subtracted  out.  This  exposure  takes  into  account  the  instrument  bias,  so  a  separate 
bias  frame  is  not  needed.  The  Sky-Pro  software  automatically  takes  a  dark  frame  of 
equal  length  for  each  light  frame  taken.  In  order  to  save  time  during  observations,  one 
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dark  frame  of  the  appropriate  length  was  taken  which  could  later  be  subtracted  from  all  of 
the  light  frames.  Each  exposure  then  needed  only  a  light  frame. 

The  other  type  of  exposure  needed  was  a  flat  field  image.  These  images  are 
typically  taken  over  a  one  to  two  second  exposure  time.  The  exact  duration  is  not  as 
important  as  all  flat  field  exposures  being  of  the  same  length.  Flat  field  images  are 
exposures  taken  of  an  illuminated  background  to  determine  the  pixel-to-pixel  sensitivity, 
or  the  quantum  efficiency  of  the  CCD  chip.6  Depending  on  when  the  observations  were 
made  on  a  particular  evening,  the  flat  field  images  could  be  taken  of  either  the 
illuminated  inside  surface  of  the  dome,  or  the  twilight  sky.  Over  the  course  of  the 
observations,  flat  fields  were  taken  using  both  techniques.  All  of  the  flat  field  images 
(typically  three  to  five)  were  averaged  to  create  a  master  flat  field  image.  This  flat  field 
correction  is  then  applied  to  the  light  frames  after  the  dark  frames  have  been  subtracted. 

Two  observations  of  each  asteroid  observed  that  evening  were  taken.  Although 
only  one  observation  was  technically  needed,  two  were  taken  mainly  for  redundancy  and 
backup  purposes.  After  approximately  one  hour,  we  returned  to  the  first  asteroid  viewed 
so  that  we  could  repeat  the  observation  cycle.  Since  some  asteroids  rose  and  set  before 
others,  we  planned  our  observations  so  that  those  that  rose  (and  therefore  set)  earliest 
were  imaged  .first  and  then  observed  the  other  asteroids  in  a  logical  sequence. 

The  observatory  log  for  the  astrometry  was  filled  out  with  the  same  information 
as  recorded  with  photometry.  Reference  Chapter  II  for  a  more  thorough  discussion  of 
that  data.  The  only  difference  was  that  no  filter  was  used  for  any  of  the  observations. 

To  go  back  for  the  second  set  of  observations  of  an  asteroid,  we  needed  to  use  the 
coordinates  that  the  telescope  was  actually  pointing  at  when  the  image  was  taken  rather 
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than  the  Guide  Star  Catalog  coordinates.  This  is  because  the  telescope  position  was 
usually  offset  slightly  from  the  Guide  Star  Catalog  coordinates  and  we  wanted  to  return 
to  the  exact  coordinates  used  previously.  Furthermore  the  Guide  Star  Catalog  gives 
coordinates  in  the  J2000  epoch  and  the  telescope  gave  the  epoch  as  the  present  time, 
J 199 8  in  this  case. 

Verifying  Capture  of  Target  Asteroid 

After  the  new  images  were  taken,  one  of  the  earlier  images  was  opened  for  a 
visual  comparison.  Although  the  asteroid  had  moved  over  this  time  interval,  sometimes 
it  was  difficult  to  notice  the  change  since  it  didn’t  move  particularly  far.  Noticing  the 
change  was  especially  difficult  in  an  image  ‘cluttered’  with  other  stars  and  objects.  Two 
example  images  of  1035  Amata  (indicated  by  the  arrow)  are  shown  in  Fig.  5.2  and  Fig. 
5.3,  which  were  taken  at  different  times.  The  Sky-Pro  software  allowed  us  to  do  a  blink 
comparison  of  the  two  frames.  A  star  common  to  both  frames  was  selected  and  its  pixel 
coordinates  (not  to  be  confused  with  celestial  coordinates)  in  each  frame  were  noted.  An 
offset  was  then  applied  to  one  picture  to  get  both  pictures  to  align.  By  then  rapidly 
blinking  the  images  back  and  forth,  we  could  visually  see  which  object  changed 
positions,  thus  identifying  it  as  our  asteroid.  The  target  asteroid  was  usually  quite  easy  to 
identify.  However,  in  the  interest  of  astronomy,  each  frame  was  carefully  inspected  for 
other  objects  that  had  changed  positions.  The  purpose  of  this  was  to  see  if  we  had  by 
chance  happened  to  catch  another  previously  unknown  asteroid.  None  of  our  frames 


indicated  such  a  discovery. 


Fig.  5.2  ST8  Image  of  1035  Amata.  UT:  05:46:58.  FOV:  6  X 10  Arcminutes 


Fig.  5.3  ST-8  Image  of  1035  Amata.  UT:  06:27:56.  FOV:  6  X 10  Arcminutes 

After  the  appropriate  data  reduction  we  decided  to  use  1035  Amata  as  the  asteroid 

for  the  orbit  determination.  Orbit  determination  would  require  observations  over  several 
weeks  and  Amata  was  favorably  positioned  for  such  observations.  In  addition,  it  had  a 
relatively  low  magnitude  of  about  16.  Furthermore,  the  published  orbital  elements  for 
Amata  were  found  on  the  Internet,  which  would  allow  us  to  verify  our  calculated 
elements  with  the  approved  elements.  In  total,  six  nights  of  observations  were  collected 
on  1035  Amata.  In  addition  to  the  author,  data  was  also  collected  by  the  Physics  480 
class  at  the  Air  Force  Academy.  A  summary  of  the  six  observation  nights  is  presented  in 
Table  5.1.  Note  the  start  time  refers  to  the  first  image  of  Amata  and  the  stop  time  refers 
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to  the  final  observation  of  Amata,  not  necessarily  of  the  entire  evening.  The  right 
ascension  and  declination  refer  to  the  telescope  pointing  coordinates,  not  Amata 


specifically. 


Date 

Start  time 

Stop  time 

#of 

Obs 

RA 

J1998 

Dec 

J1998 

UT  98  Jan  21 

UT  5:48 

UT  6:35 

B 

3  44  42.4 

42  12  32 

UT  98  Jan  22 

UT  5:12 

UT  6:44 

3  45  12.5 

42  08  35 

UT  98  Jan  27 

UT  5:10 

UT  6:46 

H 

3  45  32.8 

41  42  54 

UT  98  Jan  28 

UT  5:07 

UT  6:30 

■a 

41  38  00 

UT  98  Feb  13 

UT  1 :37 

UT  3:26 

B 

40  30  05 

UT  98  Mar  13 

UT  2:11 

UT  3:23 

6 

4  18  21.0 

39  20  08 

Table  5.1.  Astrometry  Observing  Log.  RA  given  in  hour,  min,  sec.  Dec  given  in  °, 


CHAPTER  VI 


1035  AMATA  DATA  REDUCTION 


Image  Corrections 

After  all  of  the  observations  have  been  collected,  initial  processing  is  performed. 
This  involves  the  dark  frame  and  master  flat  field  image  discussed  in  the  previous 
chapter.  These  procedures  are  performed  with  the  Sky-Pro  software.  The  dark  frame  is 
simply  subtracted  from  each  light  frame.  What  is  actually  being  subtracted  is  the  exact 
pixel  value  of  each  pixel  in  the  dark  frame  from  the  corresponding  pixels  in  the  light 
frame.  Flat  field  corrections  are  applied  in  a  more  complex  manner.  For  visualization 
purposes,  the  flat  field  can  be  thought  of  as  being  divided  into  the  light  frames.  An 
example  of  a  flat  field  image  can  be  seen  if  Fig.  6.1.  It  is  usually  good  practice  to  save 
both  the  raw  images  and  the  corrected  images  should  the  processing  need  to  be 
reaccomplished  for  whatever  reason. 


Fig.  6.1  ST-8  Flat  Field  Image 
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Astrometry  Software 

After  this  initial  processing  the  images  are  ready  to  be  processed  to  determine  the 
right  ascension  and  declination.  The  software  used  to  do  this  is  called  Astrometrica, 
version  3.1b.12  Astrometrica  is  a  shareware  program  but  the  author  requests  purchase  for 
extended  use  and  full  capabilities.13  The  remainder  of  this  chapter  will  summarize  the 
operation  of  Astrometrica,  its  capabilities,  credentials,  and  the  mathematical  techniques  it 
employs.  It  should  be  noted  that  for  our  purposes  Astrometrica  was  used  to  perform 
astrometry  of  asteroids,  but  it  can  also  be  used  for  comets,  moons,  and  variable  stars  to 

name  a  few  of  its  other  uses.  For  more  information,  visit  either  the  Astrometrica  home 

1  ^  •  1 
page  or  the  Astrometrica  user  manual. 

Astrometrica  was  originally  created  using  Borland  Pascal  7.0®.  It  was  created  to 
be  a  useful  tool  for  both  amateur  and  professional  astronomers  to  perform  CCD 
astrometry.  It  is  currently  used  in  30  countries  on  seven  continents  and  has  been  used  for 
such  purposes  as  discovering  a  new  satellite  of  Uranus  and  the  first  Allen-type  asteroid 
discovered  by  an  amateur  astronomer.  The  Astrometrica  home  page  identifies  the  Winter 
1995  issue  of  “CCD  Astronomy”  magazine,  pp.  20-22,  as  being  a  source  of  further 
information.12 

When  first  using  Astrometrica  the  location  of  the  observatory  must  be 
programmed.  This  data  is  used  to  calculate  the  coordinates  of  the  asteroid  and  is  included 
in  the  report  generated  by  Astrometrica.  The  corrected  ST-8  images  are  first  opened  into 
Astrometrica.  Astrometrica  is  capable  of  reading  several  different  data  formats, 
including  Santa  Barbara  Instrument  Group  (SBIG),  the  format  used  by  the  ST-8.  If 
necessary,  a  smoothing  or  median  filter  function  can  be  used  to  correct  any  galactic  or 
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cosmic  ray  pixel  spikes.  This  function  is  useful  for  calculating  magnitudes,  but  since  that 
was  not  the  primary  objective  of  our  observations,  it  was  not  performed  on  any  of  our 
images. 

Performing  Astrometry 

The  first  step  in  the  data  reduction  is  to  identify  the  background  brightness.  A 
box  appears  on  the  screen,  which  is  moved  over  an  area  with  no  stars  as  close  to  the 
target  object  as  possible.  Next,  the  position  and  brightness  of  the  target  object  is 
identified  using  two  similar  boxes.  The  box  size  should  be  adjusted  to  encase  the  entire 
object,  yet  being  as  close  a  fit  as  possible.  The  purpose  for  identifying  the  background 
brightness  is  similar  to  that  for  obtaining  a  dark  frame  exposure.  The  average  brightness 
value  inside  the  box  is  subtracted  from  the  target  object  and  reference  stars. 

After  identifying  the  target  object,  the  reference  or  comparison  stars  need  to  be 
identified.  Astrometrica  uses  the  Hubble  Guide  Star  Catalog  as  its  source  of  reference 
stars  which  is  convenient  for  our  purposes  since  we  used  that  software  in  locating  the 
asteroid  for  observations.  A  library  image  is  opened  on  the  screen,  which  identifies  the 
stars  for  which  Astrometrica  (via  the  Guide  Star  Catalog)  knows  the  locations.  It  is 
necessary  to  determine  which  stars  in  the  library  image  correspond  to  the  stars  in  our 
CCD  images.  A  box  similar  to  the  one  used  earlier  is  used  to  select  the  reference  stars. 
At  least  four  are  highly  recommended  and  seven  to  twelve  are  preferable,  provided  there 
are  that  many  in  all  of  the  CCD  images  and  that  they  do  not  have  close  companion  stars. 
It  is  necessary  to  identify  the  reference  stars  in  both  the  library  image  and  the  CCD 
image,  being  careful  to  select  them  in  the  same  order  each  time. 
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A  centroid  for  the  object  and  each  reference  star  is  calculated  based  on  a  center  of 
gravity  calculation  using  the  number  of  photons  collected  on  the  CCD  pixels.  A  contrast 
index  is  then  calculated  from  a  linear  regression.  This  index  is  used  to  create  a  brightness 
ratio  for  each  reference  star,  which  is  compared  against  the  Guide  catalog  information. 
The  data  from  each  reference  star  is  compiled  into  one  single,  virtual  reference  star, 
which  is  used  to  calculate  the  magnitude  of  the  target  object.  The  coordinates  are 
calculated  by  means  of  a  least  squares  fit.  This  produces  a  set  of  plate  constants  that  can 
be  translated  into  rectangular  coordinates.  A  series  of  coordinate  transformations  then 
adjusts  these  coordinates  to  right  ascension  and  declination. 

The  only  other  user  inputs  required  are  some  of  the  observing  log  parameters  such 
as  asteroid  name,  exposure  duration,  and  the  universal  time  at  the  midpoint  of  the 
exposure.  Astrometrica  calculates  the  right  ascension,  declination,  and  magnitude  from 
this  information.  The  precision  of  the  right  ascension  given  is  0.001  s  and  the  precision 
of  the  declination  is  0.01”,  although  the  final  digit  is  not  significant  for  absolute 
positions.  The  results  are  outputted  using  a  report  file  in  a  format  preferred  by  the  Minor 
Planet  Center. 

Reporting  Observations  to  the  Minor  Planet  Center 

By  adding  appropriate  header  information,  this  report  can  be  e-mailed  to  the 
Minor  Planet  Center,  which  then  adds  the  observations  to  their  database  to  update  the 
ephemeris  on  that  asteroid.  For  more  information  on  the  report  format,  reference  the 
Minor  Planet  Center  Internet  site.14  Each  observatory  is  given  an  identification  number 
that  is  published  in  the  Minor  Planet  Center  circulars  and  Internet  web  pages  to 
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acknowledge  the  contributions  from  that  observatory.  The  code  for  the  U.S.  Air  Force 
Academy  Observatory  is  712  and  can  be  seen  on  the  Minor  Planet  Center  Internet  site.10 
All  of  our  observations  were  submitted  in  two  separate  data  batches  to  the  Minor  Planet 
Center.  An  example  of  one  of  these  submissions  can  be  found  in  Appendix  K.  It  is  not 
in  the  specific  format  requested  by  the  Minor  Planet  Center,  but  it  is  summarized  to  make 
it  easier  to  understand. 

The  validity  of  all  observations  sent  to  the  Minor  Planet  Center  is  confirmed  by 
comparing  them  to  the  existing  data  for  that  asteroid.  Gareth  Williams,  the  Associate 
Director  of  the  Minor  Planet  Center,  explained  in  an  e-mail  that  all  of  the  observations 
are  propagated  at  a  nearby  200-day  epoch.15  The  resulting  orbit  is  checked  to  see  if  it 
falls  within  the  tolerances  of  their  existing  orbit  data.  The  exact  tolerance  was  not 
specified.  A  return  e-mail  was  sent  to  us  indicating  which  of  our  observations  had 
unacceptable  error.  Of  all  the  observations  submitted  by  the  USAFA  Observatory,  only 
one  was  returned  as  invalid  and  that  was  for  a  different  asteroid  than  Amata.  The  reason 
for  the  error  in  that  observation  could  possibly  be  do  to  slight  errors  in  the  Guide  Star 
Catalog  coordinates  for  the  comparison  stars,  but  the  exact  reason  is  unknown.  The 
confirmation  of  the  validity  of  our  data  is  significant.  It  verifies  that  the  observational 
and  data  reduction  procedures  followed  were  accurate  and  that  the  object  we  suspected 
was  1035  Amata  was  in  fact  that  object.  Verifying  the  validity  of  the  data  allowed  us  to 
proceed  in  determining  the  orbital  elements. 


CHAPTER  VII 


1035  AMATA  ORBIT  DETERMINATION 

Gauss’  Method  of  Angles  Only  Orbit  Determination 

Thirty-two  observations  of  1035  Amata  were  collected  from  six  nights  of 
observations.  These  data  are  presented  in  Appendix  C.  The  data  included  are  the  year; 
month;  day  and  decimal  part;  right  ascension  in  hours,  minutes,  seconds;  and  declination 
in  degrees,  minutes,  seconds.  The  computed  absolute  magnitudes  and  the  filter  used  for 
the  observations  are  also  included.  Since  none  of  the  astrometry  observations  used  a 
filter,  the  V  represents  visible.  The  magnitudes  and  filter  were  deleted,  however,  before 
the  file  was  used  by  the  orbit  determination  software.  The  numbers  preceding  certain 
lines  indicate  which  input  files  those  data  were  used  in.  For  example,  the  three  lines  with 
a  preceding  T  were  all  used  for  the  first  input  test  case  in  the  order  of  earliest  to  latest 
observation  date.  These  numbers  are  for  reference  only  and  were  not  used  in  the  actual 
input  file. 

Advising  in  the  orbit  determination  for  1035  Amata  was  Mr.  Roger  Mansfield. 
He  recommended  that  Gauss’  method  of  angles  only  orbit  determination  would  be 
appropriate  for  this  problem  since  the  data  included  three  position  measurements  (angles) 
and  the  times  of  the  observations.  This  method  had  an  equation  with  six  unknowns, 
which  must  be  solved  to  determine  the  orbit.  Since  each  observation  provides  two 
angles,  a  minimum  of  three  observations  are  required.  Unless  otherwise  noted,  the 
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following  summary  of  Gauss’  angles  only  method  of  orbit  determination  is  taken  from 
J.M.A.  Danby’s  Fundamentals  of  Celestial  Mechanics .'6 

The  first  consideration  in  using  this  approach  is  that  for  some  problems  a  solution 
cannot  be  found.  This  can  result  from  inaccurate  position  measurements  or  an 
unfavorable  geometry  in  the  relative  positions  of  the  observations  to  each  other. 
Additionally,  for  some  deep  space  objects  with  peculiar  orbits,  an  orbit  may  simply  not 
be  determinable  with  this  method.  Gauss  assumes  that  the  three  observations  lie  in  a 
plane.  This  assumption  would  not  be  valid  if  third-bodies  were  exerting  perturbing  forces 
on  the  asteroid,  which  would  change  its  orbit  plane.  Although  several  thousand  asteroids 
have  orbits  close  to  Amata  (celestially  speaking),  the  small  mass  of  these  objects  is  not 
likely  to  exert  a  large  perturbing  force.  Perturbations  from  larger  objects  such  as  Jupiter 
or  Mars  are  probably  not  a  concern  over  the  short  time  interval  (three  months)  during 
which  the  observations  were  made.  The  relative  positions  of  Earth,  Jupiter,  Mars  and 
Amata  at  the  time  of  observation  are  unknown. 

Since  all  the  vectors  lie  in  the  same  plane,  it  is  possible  to  define  one  of  them,  the 
middle  observation,  in  terms  of  the  other  vectors  with  two  scalar  constants. 


r2  =  cir,  +  c3r3  (7.1) 

One  of  the  most  difficult  aspects  of  Gauss’  method  is  selecting  the  initial  values  for  these 
scalars.  Danby  recommends  using  the  opening  terms  from  the  f  and  g  series  to  select 


these  values. 
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Cl  =  g 3  /  (flg3  -  glfj) 

C3  =  -  gl  /  (flg3  -  glf3) 

After  successive  iterations,  the  constants  are  updated,  but  an  accurate  initial  guess  is 
crucial  in  determining  whether  or  not  the  process  will  converge. 

These  constants  are  then  used  in  calculating  Gauss’  sector-triangle  ratios.  This  is 
a  ratio  of  the  entire  area  swept  out  by  the  asteroid  orbit  between  observations  to  the  area 
of  the  triangle  swept  out  between  observations.  The  concept  behind  this  calculation  is 
Kepler’s  Second  Law,  which  equates  equal  areas  in  equal  times.17  An  illustration  of  this 
ratio  is  provided  in  Fig  7.1.  Sector-triangle  ratios  are  computed  for  observations  1  to  2,  2 
to  3,  and  1  to  3.  Fig.7.1  shows  the  calculation  for  observations  1  to  3.  A  series  of 
coordinate  transformations  are  then  performed  to  convert  the  observations  from  the 
equatorial  system  to  a  new  system.  A  similar  transformation  is  performed  on  the  Sun 
coordinates  to  be  in  the  same  reference  system  as  the  observations.  After  this,  the 
geocentric  distances  to  the  asteroid  are  calculated  as  follows 

P2  =  (-CiZj  +  Z2  -  c3Z3)  /  o2  (7 .4) 

p3  =  (p2ft2  +  CiYi  -  Y2  +  c3Y3)  /  c3p3  (7.5) 

Pi  =  (P2^-2  -  C3P3A-3  +  C]Xi  -  X2  +  C3X3)  /  Cl  (7.6) 

where  X,  p,  and  o  are  components  of  the  unit  vectors  in  the 
direction  of  the  observations  and  X,  Y,  and  Z  are  the  solar 


(7.2) 

(7.3) 


coordinates. 


41 


The  process  is  repeated  until  the  change  in  the  new  geocentric  distances  is  less  than  a 
specified  tolerance.  From  the  final  range  distance,  the  desired  orbital  elements  are 
calculated. 

r3 


rl 

ju||  Sec-Tri  Ratio  =  (Area  1  +  Area  2)  /  Area  2 

Fig,  7.1  Example  Sector  Triangle  Ratio 

GEM  Software 

It  was  decided  to  use  the  set  of  computer  programs  for  calculating  orbital 
elements  that  came  with  Danby’s  Fundamentals  of  Celestial  Mechanics .16  The  specific 
program  of  interest  was  the  Gauss-Encke-Merton  (GEM)  method,  named  after  its  three 
main  contributors.  The  Turbo  Pascal  version  of  this  program  was  used.  To  begin,  the 
unit  vectors  for  the  direction  cosines  of  the  observations,  and  the  Sun  coordinates  at  the 
observation  times  are  calculated.  The  exact  Sun  coordinates  are  calculated  based  on  the 
Julian  Date  of  the  observation  and  tabular  information  provided  by  Bretagnon  and  Simon. 
These  are  both  calculated  in  the  ecliptic  plane  with  an  obliquity  of  the  ecliptic  unique  for 
each  observation  Julian  Date. 

By  Danby’s  admission,  the  GEM  program  is  not  a  completely  accurate  program 
since  it  does  not  account  for  planetary  aberration.  Danby  also  indicates  that  more 
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accurate  Sun  position  angles  should  be  determined.18  Planetary  aberration  indicates  that 
the  time  an  observation  is  made  does  not  indicate  the  true  position  of  an  object  for  that 
recorded  time  since  the  light  reflected  from  it  travels  to  the  Earth  with  a  finite  speed. 
Despite  its  potential  shortcomings,  the  GEM  program  was  attempted  and  appropriate 
modifications  and  corrections  were  made  where  needed. 

Some  of  the  initial  modifications  made  were  to  add  two  additional  data  input 
procedures,  which  allowed  the  data  to  be  entered  either  from  the  keyboard  as  the  program 
was  running  or  to  read  data  from  an  external  input  file.  Additionally  a  declination 
conversion  was  added  to  allow  the  declination  to  be  input  as  degrees,  minutes,  seconds 
and  then  converted  to  degrees  and  fractional  degrees.  This  is  the  format  provided  by  the 
astrometry  process  and  most  commonly  used  by  astronomers.  A  user  query  was  also 
added  after  each  iteration  loop  to  ask  the  user  if  it  should  continue  with  the  next  iteration. 
The  data  printed  to  the  screen  at  each  query  allowed  to  user  to  see  if  the  program  was 
converging  or  not.  The  specific  values  checked  for  convergence  were  the  geocentric 
distances  and  the  sector-triangle  area  ratios.  There  were  three,  separate  values  for  each  of 
these,  one  for  each  input  observation.  All  external  procedures  were  also  incorporated 
into  the  main  GEM  program  file  so  separate  procedure  files  were  not  needed.  Comments 
and  structure  changes  were  also  made  throughout  to  make  the  program  easier  to  read  and 
understand.  The  GEM  program  with  all  of  our  modifications  can  be  found  in  Appendix 
D. 

Since  the  data  included  thirty-two  observations  of  Amata,  six  sets  of  random 
combinations  of  three  observations  each  were  created.  In  order  to  ensure  sufficient  time 
and  distance  between  observations,  multiple  observations  from  the  same  evening  were 
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not  used.  In  general,  an  observation  was  selected  from  the  beginning,  middle,  and  end  of 
the  data  list.  The  actual  data  sets  used  are  numbered  in  Appendix  C. 

The  program  did  converge  for  these  data  and  was  close,  but  not  identical  to  the 
orbital  elements  published  by  the  Minor  Planet  Center.  The  published  orbital  elements 
from  the  Minor  Planet  Center  can  be  found  in  Table  7.1.  Mr.  Gareth  Williams  explained 
that  the  Minor  Planet  Center  uses  several  different  orbit  determination  algorithms,  but 
most  of  them  are  a  derivative  of  Gauss’  method.15  The  calculated  elements  from  six  test 
cases  using  the  GEM  program  can  be  found  in  Table  7.2.  The  complete  output  files  from 
the  GEM  program  can  be  found  in  Appendix  E.  To  find  errors  in  the  program,  hand 
calculations  were  performed  using  the  Amata  data.  An  error  in  the  calculation  of  one  of 
the  arrays  was  found  in  the  program.  Correcting  it  did  not  significantly  improve  the 
results  however,  since  the  error  was  in  the  direction  of  the  K  unit  vector  and  the 
inclination  in  this  direction  is  less  than  20°.  A  modified  subroutine  to  calculate  the  Sun 
angles  that  converted  all  the  elements  to  the  J2000  epoch  was  also  tried  but  it  failed  to 
significantly  improve  the  results.  No  attempt  to  account  for  planetary  aberration  was 
made  since  it  was  felt  that  this  too  would  provide  only  minimal  improvements  in  the 
accuracy. 


Minor  Planet  Center 


3.137178  :  1&08732 

P  (years)  5-56  M  (°)  85.82541 

IllpS#  0.2026701  2.20159 

0.17737536  l^^jpaB  323.12242 
Table  7.1  Minor  Planet  Center  Orbital  Elements  19 
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I  GEM 

Input  1 

Input3|S& 

Input  3  Input  4 

Inputs 

Input  6 

a  (AD) 

3.2182661 

3.2101537 

3.2268884  3.2727042 

3.2071182 

3.2503599 

P  (years) 

5.7735259 

5.7517096 

5.7967439  g  5.9206355 

5.7435533 

5.8601047 

0.1597266 

0.1632919 

0.1559225  0.1237473 

0.1660712  f 

0.1431401 

n  (°/day) 

0.1707148 

0.1713623 

0.1700310  i  0.1664731 

0.1716057  | 

0.1681926 

Inch  (°) 

li  m-  ss\\  7:;;' 

18.062012 

18.063142 

18.061130  18.049346 

18.061310  | 

18.053590 

M(°) 

Q(°) 

1.4028082 

1.4833787 

fi:  ,  '■  , 

1.3648144  0.9540862 

1.5977124 

1.1920249 

0)  (°) 

329.27072 

328.62022 

329.91056  334.45223 

328.51103 

332.18644 

Table  7.2  GEM  Orbital  Elements.  J2000  Epoch 

ORBDET  Software 

It  was  then  decided  to  try  the  ORBDET  orbit  determination  program  by 
Montenbruck  and  Pfleger.17  Without  any  modifications,  this  program  yielded  results  that 
were  significantly  closer  to  the  published  values  than  the  GEM  program  did.  The  orbital 
elements  calculated  by  ORBDET  from  the  six  test  cases  can  be  found  in  Table  7.3.  A 
slight  modification  to  the  input  data  needed  to  be  made  for  the  ORBDET  program.  The 
decimal  part  of  the  day  was  multiplied  by  24  to  give  hours  and  decimal  hours.  An 
example  of  an  input  file  for  ORBDET  can  be  found  in  Appendix  G.  The  complete  output 
files  from  ORBDET  can  be  found  in  Appendix  H.  ORBDET  also  accounted  for 
planetary  aberration. 
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ORBDET 

Input! 

•  LOipilt  2 

Input  3 

Input  4 

Inputs  :  ? 

Input 

f.  ■  •'7  5  1 

a  (AU) 

3.139193 

3.139943 

-0.733254  , 

3.140728 

3.140653 

3.141972 

P  (years) 

5.5621 

::|i|!5.5640 

0.6279  | 

5.5661 

5.5659 

1 

5.5694 

0.202770 

0.202402 

4.394909  | 

0.202682 

a. 4jr  ; 

0.202028 

0.201908 

n  f/day) 

0.177205 

jjgg| 

1.569720  ' 

0.177076 

0.177082  j 

■"  'V  ;  :  :A 

0.176970 

Inch  (°) 

18.0870 

18.0866 

18.6437 

18.0868 

18.0864 

P 

18.0861 

M(°) 

85.670976 

85.671276 

-332.9817  f 

85.575542 

85.670941 

85.595332 

■  ■  ■■ 

Q(°) 

2.1808  ! 

2.1676 

340.4056 

2.1664 

2.1577  2.1416 

o>  (°) 

323.2704  | 

323.3218 

158.9076  I 

323.3803 

323.3711  1 

232.4682 

Table  7.3  ORBDET  Orbital  Elements.  J2000  Epoch 

With  this  information,  we  decided  to  abandon  further  attempts  at  debugging  the 

GEM  program.  In  general,  the  mathematical  techniques  employed  in  Danby’s  text  are 
thought  to  be  more  robust  than  those  used  by  Montenbruck  and  Pfleger.  Montenbruck 
and  Pfleger  however,  are  more  concerned  with  reliable  software  development.  Therefore 
even  though  mathematically  not  as  robust  as  GEM,  ORBDET  has  been  properly 
debugged,  yielding  results  that  are  more  accurate. 

The  following  description  of  the  ORBDET  program  comes  from  the  Montenbruck 
and  Pfleger  text.17  ORBDET  uses  Bucerius’  method,  which  is  derived  from  Gauss’ 
method  of  angles  only  orbit  determination  but  with  more  simplifications.  The  heart  of 
Gauss’  method,  the  sector-triangle  ratio  calculations,  is  used  here  as  well  as  in  GEM. 
The  main  simplification  in  Bucerius'  method  to  Gauss’  method  is  that  it  is  unable  to 
account  for  orbital  paths  crossing  a  tear-drop  shaped  region  of  the  space  around  the  Sun 
and  inside  the  Earth’s  orbit  defined  by  Charlier’s  boundary  line.  fig.  7.2  illustrates  this 
region.  Observations  from  this  region  could  have  two  possible  solutions  which  Bucerius' 
method  cannot  handle.  An  object  with  an  orbit  crossing  this  region  might  be  a  near  earth 
asteroid  or  comet  with  a  high  eccentricity.  However,  since  the  orbit  of  Amata  never  takes 
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it  inside  this  region,  it  can  always  be  observed  at  opposition,  which  the  Bucerius’  method 
is  capable  of  handling.  Although  ORBDET  cannot  handle  these  situations,  it  does  flag 
when  such  a  situation  occurs.  This  would  allow  the  user  to  attempt  a  more  complete 
method  of  orbit  determination. 


Fig.  7.2  Charlicr’s  Boundary  17 

Additional  calculations  for  other  orbital  elements  were  added  to  ORBDET  to 
allow  us  to  compare  with  the  Minor  Planet  Center  elements.  Those  equations  for  the 
mean  motion,  mean  anomaly,  and  the  period  are 


n  =  kg  *  a  <-3/2) 

(7.7) 

M  =  n  *  (JD2  —  JDmpc) 

(7.8) 

P  =  (360  /  n)  *  (1  /  365.25) 

(7.9) 

where  kg  is  the  Gaussian  constant,  JD2  is  the  JD  for  the  second 
observation,  and  JDmpc  is  the  Minor  Planet  Center  epoch 

Notice  that  the  mean  motion  was  corrected  to  the  epoch  used  by  the  Minor  Planet  Center 
to  allow  for  a  direct  comparison  of  this  parameter.  The  main  program  that  calculates  the 
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orbital  elements  can  be  found  in  Appendix  F.  ORBDET  uses  a  large  number  of  functions 
and  procedures  found  in  external  files.  To  see  these  files  reference  the  software 
associated  with  Montenbruck  and  Pfleger’s  text. 17 

FINDORB  Software 

A  last  minute  discovery  was  made  in  the  preparation  of  this  thesis.  A  batch  least 
squares  program  was  found  on  the  Internet,  which  was  capable  of  handling  an  entire  list 
of  observations.  The  name  of  this  program  is  FIND  ORB,  written  by  Bill  Gray  and 
maintained  at  the  Project  Pluto  web  page.20  FIND_ORB  is  ideal  software  for 
astronomers  who  submit  data  to  the  Minor  Planet  Center.  The  output  format  is  identical 
to  that  used  in  the  Minor  Planet  Center  circulars.  The  epoch  used  by  the  Minor  Planet 
Center  can  also  be  directly  input  to  allow  direct  comparison  of  all  the  computed  elements. 

FIND  ORB  initially  assumes  a  distance  of  1  AU  from  the  Sun  to  the  first  and  last 
observations.  Through  iterations,  this  distance  is  modified.  The  program  will  converge 
faster  if  a  more  accurate  initial  distance  guess  is  input.  This  method  is  referred  to  as  the 
Herget  method,  which  is  used  to  get  an  initial  estimate  of  the  orbit.  A  residual  is 
calculated  for  each  iteration,  which  is  the  observation  values  minus  the  calculated  values. 
The  Herget  step  is  repeated  until  the  residuals  aren’t  decreasing  anymore.  With  this 
estimation  of  the  orbit  a  full  step  is  taken  which  involves  the  least  squares  “best  fit” 
process.  A  few  full  steps  are  performed  until  the  residuals  fail  to  decrease  any  further. 

FIND  ORB  also  has  the  capability  of  allowing  for  perturbations  of  all  the  planets 
and  the  moon,  although  each  body  selected  increases  the  computing  time  required  for 
each  iteration.  It  is  interesting  that  including  all  of  the  possible  perturbations  did  not 
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decrease  the  residuals  any  further.  This  would  indicate  that  the  two-body  assumption  for 
Gauss’  method  was  valid.  A  possible  reason  third-body  affects  weren’t  noticed  is  that  the 
observations  were  collected  over  only  a  couple  of  months.  Had  observations  been  made 
over  several  years,  perturbation  effects  probably  would  have  been  noticed. 

Bill  Gray,  the  author  of  FIND_ORB  indicates  that  residuals  of  less  than  one 
arcsecond  would  indicate  quite  accurate  observations.  For  our  observations, 
FIND  ORB  calculated  a  final  residual  of  0.360  arcseconds,  which  is  another 
confirmation  of  the  validity  of  our  observational  and  data  reduction  procedures.  The 
calculated  elements  matched  the  MPC  values  to  about  four  or  five  decimal  places  for 
most  parameters. 

The  average  values  of  the  computed  elements  from  each  method  used  are  shown 
in  Table  7.4.  Note  that  the  ORBDET  elements  from  test  case  3  were  not  used  because 
ORBDET  failed  to  find  a  reasonable  solution  for  the  orbit  for  that  test  case.  The  three 
programs  used  are  also  arranged  in  order  of  increasing  accuracy. 
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Table  7.4  Average  Orbital  Elements 

CHAPTER  VIII 


CONCLUSIONS  AND  SUGGESTIONS  FOR  FURTHER  RESEARCH 

Photometry  Results 

Differential  photometry  proved  to  be  an  effective  method  for  calculating  the 
rotational  period  of  an  asteroid.  To  restate  the  results  of  our  work,  a  9.210  ±  0.005-hour 
synodic  period  was  calculated  for  583  Klotilde.  We  have  a  high  level  of  confidence  in 
our  method  and  procedures  at  arriving  at  this  period.  Of  course  with  any  observational 
procedure,  there  is  the  possibility  of  error  especially  considering  that  visual  judgments 
were  made  for  some  of  our  conclusions.  We  do  not  feel  however,  that  these  judgments 
have  invalidated  our  conclusions. 

Additional  Photometry  Research 

There  are  still  numerous  asteroids  in  our  solar  system  in  a  variety  of  orbits  both 
close  to  and  distant  from  the  Earth  that  have  not  had  their  rotational  periods  calculated. 
Additional  asteroids  have  had  preliminary  work  done,  but  no  final  rotational  period  has 
been  determination.  The  procedures  outlined  in  this  thesis  could  be  duplicated  for  many 
of  those  asteroids. 

Of  particular  interest  in  our  research  was  the  indication  of  discovering  a 
previously  unknown  variable  star.  Further  observations  and  analysis  are  needed  to  verify 
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our  initial  data  and  characterize  the  nature  of  this  star’s  variability.  The  star  in  question  is 
GSC  223:1761  (Guide  Star  Catalog  designation)  located  in  the  constellation  Cancer. 

Astrometry  Results 

The  observational  techniques  and  data  reduction  procedures  used  to  determine 
asteroid  positions  were  verified  by  confirmation  from  the  Minor  Planet  Center  as  well  as 
the  ORBDET  and  FINDORB  orbit  determination  programs.  The  results  obtained 
through  ORBDET  closely  matched  the  elements  calculated  by  the  Minor  Planet  Center 
and  the  batch  least  squares  method  used  by  FIND_ORB  yielded  even  closer  results.  The 
GEM  program  also  has  potential  to  be  an  accurate  method  of  orbit  determination,  but 
further  debugging  needs  to  be  done  on  it  first.  Although  Gauss’  method  of  angles  only 
orbit  determination  worked  for  our  data  there  may  be  orbits  where  this  method  will  not 
converge  and  it  will  not  be  possible  to  determine  an  object’s  orbit. 

Additional  Astrometry  and  Orbit  Determination  Research 

There  are  also  several  areas  for  further  research  in  astrometry.  The  Minor  Planet 
Center  continuously  updates  its  critical  list  of  asteroids.  They  are  always  interested  in 
observations  from  amateur  astronomers  to  augment  their  ephemerides.  Coordinates  can 
be  computed  for  any  of  these  asteroids  and  reported  to  the  Minor  Planet  Center  without 
having  to  perform  the  orbit  determination  done  in  the  project.  Similar  work  can  also  be 
accomplished  for  comets.  There  are  undoubtedly  hundreds  of  other  minor  planets  and 
comets  in  our  solar  system  that  have  not  yet  been  discovered.  An  observant  astronomer 
might  be  fortunate  enough  to  discover  a  new  comet  or  asteroid. 
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The  orbit  determination  methods  used  here  also  suggest  areas  of  additional  work. 
ORBDET  was  only  modified  to  calculate  additional  orbital  elements  for  comparison  with 
the  Minor  Planet  Center  elements.  Since  our  database  of  observations  was  quite  large, 
the  program  could  be  modified  to  read  in  several  observations  and  apply  a  batch  least 
squares  approach  to  average  them.  A  similar  process  could  be  accomplished  with  GEM 
once  all  the  bugs  in  the  program  have  been  removed. 
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Differential  Photometry  Spreadsheets 
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APPENDIX  B 

Submitted  Rotational  Period  Paper 


CCD  PHOTOMETRY  OF  583  KLOTILDE  AT  THE 
U.S.  AIR  FORCE  ACADEMY  OBSERVATORY 


Daniel  C.  Burtz 

Master  of  Engineering  Program  Office 
University  of  Colorado  at  Colorado  Springs 
Colorado  Springs,  CO  80918 

Charles  J.  Wetterer 
Department  of  Physics 
United  States  Air  Force  Academy 
USAF  Academy,  CO  80840 

CCD  photometry  observations  of  asteroid  583  Klotilde,  taken 
during  1998  February  and  March  at  the  U.S.  Air  Force  Academy 
Observatory  are  reported.  The  asteroid  was  observed  as  part  of  a  masters 
degree  thesis  project.  A  synodic  period  of  9.210  ±  0.005  hours  was 
determined  from  four  nights  of  observations.  The  observed  amplitude  was 
0.18  magnitudes. 


All  observations  were  made  at  the  U.  S.  Air  Force  Academy  Observatory.  A 
Photometries  (PM512)  CCD  camera  attached  to  a  61 -cm  Cassegrain  telescope  was  used 
to  take  five  minute  exposures  of  the  asteroid  through  a  standard  Johnson  R-band  filter. 
All  images  were  processed  using  NOAO’s  IRAF  package  and  differential  photometry  of 
the  asteroid  with  respect  to  nearby  stars  was  performed.  Differential  photometry 
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measurements  between  several  stars  in  each  field  were  conducted  to  ensure  the 
comparison  stars  were  not  variable.  The  data  revealed  that  one  of  the  comparison  stars 
from  one  night  was  in  fact  variable.  This  star  was  not  used  as  a  comparison  star.  Since  it 
is  not  identified  as  a  variable  star,  further  observations  are  currently  being  conducted  to 
determine  the  nature  of  its  variability. 

Potential  asteroid  targets  were  selected  by  first  setting  an  upper  magnitude  limit 
for  asteroids  at  a  favorable  viewing  position  during  February  and  March.  Information 
provided  by  Harris  (1997)  was  then  consulted  to  see  which  of  these  asteroids  did  not  have 
previously  determined  rotational  periods.  The  asteroid  583  Klotilde  was  chosen  as  our 
primary  target. 

Asteroid  583  Klotilde  was  observed  48  times  over  a  5.5  hour  period  (4.8  to  10.3 
hours  UT)  on  UT  98  Feb  13;  52  times  over  a  5.6  hour  period  (3.4  to  9.0  hours  UT)  on  UT 
98  Feb  22;  56  times  over  a  7.5  hour  period  (1.9  to  9.4  hours  UT)  on  UT  98  Mar  12;  and 
50  times  over  a  4.7  hour  period  (2.8  to  7.5  hours  UT)  on  UT  98  Mar  13.  The  data  from 
UT  98  Mar  12  suggests  a  triple  minima  lightcurve  with  each  minimum  having  a 
distinctive  depth.  Combining  this  night  with  the  next  and  visually  inspecting  the 
resulting  lightcurve  indicates  a  synodic  period  of  9.2  ±0.1  hours  containing  three 
minima.  Adding  the  observations  from  the  other  nights  refines  the  synodic  period  to 
9.210  ±  0.005  hours.  The  composite  lightcurve  using  this  period  is  shown  in  Figure  1. 
The  observed  amplitude  was  0.18  magnitudes. 
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Figure  Captions:  (For  figure  submitted  see  Fig.4.5.) 

Figure  1.  Composite  lightcurve  for  583  Klotilde  based  on  a  9.210  hour  synodic  period. 
Zero  phase  =  UT  98  Mar  13,  00:00:00.  Magnitudes  calibrated  roughly  with  the  Hubble 
Guide  Star  Catalog  magnitudes  (Russell  it  al.  1990)  of  comparison  stars,  adjusted  to 
common  Earth  and  Sun  distances  (1  A.U.  for  both),  and  scaled  to  UT  98  Mar  13  data. 
Lightcurve  is  corrected  for  light  travel  time. 


APPENDIX  C 


Computed  Amata  Coordinates  /  Input  Data 


Test  case,  Year,  Month,  Day.day,  RA  RA  RA.ra,  Dec  Dec  Dec.dec,  magnitude,  filter 

1  1998  01  21.24164  03  44  50.43  42  12  41.6  15.9  V 

1998  01  21.24550  03  44  50.41  42  12  40.7  15.9  V 

1998  01  21.27009  03  44  50.55  42  12  32.5  16.0  V 

1998  01  21.27427  03  44  50.54  42  12  31.2  16.0  V 

2  1998  01  22.21771  03  44  55.41  42  07  42.6  15.9  V 

1998  01  22.22047  03  44  55.42  42  07  41 .2  1 5.9  V 

1998  01  22.27830  03  44  55.73  42  07  23.9  16.1  V 

1998  01  22.28038  03  44  55.74  42  07  23.1  16.1  V 

3  1998  01  27.21590  03  45  46.10  41  43  02.5  16.5  V 

1998  01  27.21848  03  45  46.10  41  43  01.6  16.6  V 

1998  01  27.25006  03  45  46.54  41  42  53.6  15.7  V 

1998  01  27.25433  03  45  46.62  41  42  51.9  16.0  V 

1998  01  27.27655  03  45  46.94  41  42  45.1  16.0  V 

4  1998  01  27.28177  03  45  47.03  41  42  43.5  16.0  V 

5  1998  01  28.21384  03  46  01.17  41  38  18.7  16.1  V 

1998  01  28.21614  03  46  01.20  41  38  18.3  15.9  V 

1998  01  28.24293  03  46  01.63  41  38  10.2  16.0  V 

1998  01  28.24633  03  46  01.68  41  38  09.1  16.1  V 

1998  01  28.26623  03  46  01.96  41  38  03.8  16.0  V 

6  1998  01  28.27069  03  46  02.02  41  38  02.5  16.0  V 

1  1998  02  13.06861  03  53  21.41  40  32  49.2  16.6  V 

2  1998  02  13.07160  03  53  21.53  40  32  48.7  16.5  V 

3  1 998  02  1 3.09028  03  53  22.22  40  32  44.8  1 6.4  V 

4  1998  02  13.09264  03  53  22.33  40  32  44.0  16.4  V 

5  1998  02  13.12376  03  53  23.49  40  32  37.6  16.4  V 

6  1998  02  13.12638  03  53  23.59  40  32  37.1  16.4  V 

1  1998  03  13.09222  04  18  34.54  39  20  22.9  16.2  V 

2  1998  03  13.09472  04  18  34.71  39  20  22.4  16.4V 

3  1998  03  13.09684  04  18  34.84  39  20  22.2  16.1  V 

4  1998  03  13.13267  04  18  37.23  39  20  18.7  16.2  V 

5  1998  03  13.13863  04  18  37.63  39  20  17.6  16.3  V 

6  1998  03  13.14065  04  18  37.76  39  20  17.3  16.1  V 


APPENDIX  D 


Modified  GEM  Program 


Program  GEM  (input,  output); 

{This  program  was  originally  produced  by  J.M.A.  Danby.  Modifications 
and  alterations  were  made  by  2Lt  Dan  Burtz  in  Mar  98  for  thesis  work  for 
a  Master  of  Engineering  degree  in  Space  Operations  from  the  University  of 
Colorado  at  Colorado  Springs.  The  data  I  used  in  this  program  was  taken 
from  the  U.S.  Air  Force  Academy  16  inch  telescope  with  the  ST8  CCD  camara. 
Images  taken  were  of  asteroid  1035,  Amata,  from  21  Jan  98  through  13  Mar  98. 

Equation  numbers  in  parenthesis  can  be  found  in  Fundamentals  of  Celestial 
Mechanics  by  J.M.A.  Danby,  2nd  Ed.,  cl992.  Published  by  Wi llmann-Bell . 

This  program  reads  in  data  from  an  input  file  in  the  following  format: 

YYYY  MM  DD.DDDDD  RA  RA  RA.RA  DE  DE  DE . D 

where  RA  and  DE  represent  right  ascension  and  declination  in  their 
respective  units.  The  program  also  contains  procedures  which  are  commented 
out  that  allow  inputing  these  values  in  real  time  from  the  keyborad,  and 
for  hardcoding  them  Into  the  program.  The  data  input  data  is  echo  checked 
to  an  output  file  along  with  the  computed  orbital  elements.  Due  to  the 
fact  that  the  iterations  in  this  program  may  not  converge,  continuing 
with  each  iteration  is  specified  by  the  user  while  running  the  program. 

The  final  orbital  elements  are  also  printed  to  the  screen. 


(Notes  from  Danby) 

Basic  program  for  determining  an  orbit  from  three  observations. 
See  section  7.3.  The  GEM  method  is  used,  where  sector- tri angle 
ratios  are  found.  Allowance  has  NOT  been  made  for  planetary 
aberration.  To  turn  this  into  a  completely  practical  program, 
this  modification  should  be  made,  and  possible  more  accurate 
solar  coordinates  used.} 


PROGRAM  VARIABLES  AND  PROCEDURES 


VARIABLES 


a 

- 

Semi-major  axis  of  orbit 

AU 

argp 

- 

Argument  of  periapsis 

rad 

cl 

- 

Scalar  component  of  R  vectors 

c3 

- 

Scalar  component  of  R  vectors 

change 

- 

Change  in  sum  of  rho's  for  each  iteration 

AU 

CONTINUE 

- 

Loop  control  variable  to  continue  with  iterations 

integer 

COUNT 

- 

Iteration  number 

integer 

d 

- 

Integer  Julian  day  for  input;  Day  of  month  for  output 

integer 

dec 

- 

Declination  of  observation 

rad 

decdeg 

- 

Degree  portion  of  input  declination 

deg 

decmi n 

- 

Minute  portion  of  input  declination 

mi  n 

decsec 

- 

Second  portion  of  input  declination 

sec 

dt 

- 

e 

- f 

- 

Eccentricity 

T 

f i lei n 

_ 

Variable  designating  the  input  data  file 

f i leout 

- 

Variable  designating  the  output  data  file 

f ullday 

- 

Julian  day  with  fractional  part 

g 

gd 

_ 

Geocentric  distance,  the  magnitude  of  each  rho 

array  of 

gk 

- 

Gaussian  constant  (sqrt  of  MU  of  sun) 

he 

- 

Heliocentric  coordinates 

array 
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hour 

i 

inc 

j 

k 

kay 

m 

MA 

min 

n 

omega 

P 

Pi 

r 

r  1 
r  2 
ra 
rh 

rhol 

rho2 

rho3 

rm 

sec 

scl 

sc2 

sx 

sy 

sz 

t 

taul 

tau3 

temp 

time 

tpc 

ut 

x 

xl 

x2 

xs 

xv 

y 

yi 

y2 

y3 

ye 

ygi 

yv 

Z 

Z  V 


-  Right  ascension  hour  portion  real 

-  Loop  control  variable  integer 

-  Inclination  of  orbit  rad 

-  Loop  control  variable  integer 

-  Loop  control  variable  integer 

-  Month  number  in  output  integer 

-  Mean  anomaly  rad 

-  Right  ascension  minute  portion  min 

-  Mean  motion  rad/day 

-  Longitude  of  the  ascending  node  through  ecliptic  rad 

-  Semi-latus  rectum  AU 

-  Arc  Cos(-l)  real 

-  Right  ascension  of  observation  rad 


-  Right  ascension  seconds  portion  sec 

-  Direction  cosines  of  Solar  coordinates  array 

-Rotated  solar  coordinates  array 


-  Time  constant  (delta  t  x  gk) 

-  Time  constant  (delta  t  x  gk) 

-  The  Julian  date  of  the  observation  time  days 

-  Decimal  part  of  full  input  Julian  Day  days 

-  Heliocentric  equitorial  coordinates  for  time  1 

-  Heliocentric  equitorial  coordinates  for  time  2 


-  Sector  triangle  ratio  for 

-  Sector  triangle  ratio  for 

-  Sector  triangle  ratio  for 

-  Year  of  observation  time; 


obs  1  -  obs  2 
obs  2  -  obs  3 
obs  1  -  obs  3 


uni tless 
unitless 
unitless 


Year  of  calculated  elements  integer 


PROCEDURES 

Coord i nates_elements 

Data 

Date_Jd 

Jd_Date 

Output 

Rotation 

Sector_tri angle- ratio 
Sun 


-  Calculates  the  orbital  elements 

-  Gets  input  times,  RA,  Dec;  converts  to  working  units 

-  Converts  a  calendar  date  to  a  Julian  date 

-  Converts  a  Julian  Date  to  a  calendar  date 

-  Writes  orbital  elements  to  screen  and  file 

-  Resolves  vectors  in  Cunningham's  reference  system 

-  Uses  Gauss'  method  to  compute  the  ratio 

-  Calculates  solor  coordinates 


FUNCTIONS 

Atan2  -  Finds  the  arc  tangent  of  two  numbers 

Obliquity  -  Calculates  the  obliquity  of  the  ecliptic  based  on  the  JD 
Q  -  Calculates  Q  function  for  Gauss'  method} 
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uses  Crt,  Dos; 

Const 

pi:  real  =  3.1415926535897932; 
gk:  real  =  0.01720209895; 


<»  > 
var 


fileout,  filein: 

text ; 

COUNT,  CONTINUE: 

integer ; 

ra,  dec,  time,  sx,  sy,  sz: 

array [1. 

.3] 

of 

real ; 

rhol,  rho2,  rho3,  gd: 

array  [  1 . 

.3] 

of 

real ; 

rh,  rm,  scl,  sc2,  he: 

array  [1 . 

.3, 

1. . 

.3]  of 

ye ,  m: 

integer; 

d,  ut,  hour,  min,  sec,  fullday: 

real ; 

deedeg,  deemin,  decsec: 

real ; 

i ,  j .  k: 

integer ; 

taul,  tau3,  cl,  c3: 

real ; 

temp,  change: 

real ; 

rl,  r2 ,  kay,  dt,  xs,  ygl,  yl,  y2, 

y3: 

real ; 

xl,  x2: 

array [1 . 

•  3] 

of 

real ; 

r ,  x ,  y ,  z ,  xv ,  yv ,  zv ,  f ,  g: 

real ; 

obi,  a,  e,  inc,  omega,  argp,  tpc, 

t: 

real ; 

n,  p,  MA: 

real ; 

Procedure  Date_Jd(var  Year,  Month:  integer;  Day,  Ut:  real; 

var  Jd :  real) ; 

var 

ryear:  real; 

Begin 

ryear  :=  Year; 

Jd  :=  367* ryear  -  Int (1 . 75* (rYear  +  Int((Month  +  9)/12)))  +  Day 
+  Int (275*Month/9)  +  1721013.5  +  Ut/24; 

end ; 

{  ~> 
Procedure  Jd_date(Jd:  real;  var  Year,  Month:  integer;  var  Day:  real); 

{Finds  the  calendar  date  from  the  Julian  day  number.  Uses  the 
procedure  given  by  Jean  Meeus  in  "Astronomical  Formulas  for 
Calculators,  pp  26,  27.} 
var 

z ,  f,  a,  b,  c,  d:  real ; 
e:  integer; 

begin 

Jd  :=  Jd  +  0.5; 
z  :=  i nt ( Jd  +  0.00001) ; 
f  :=  Jd  -  z; 

if  z  <  2299161.0  then  a  :=  z 
else 

begin 

a  :=  int ( (z  -  1867216 . 25) /36524 . 25) ; 
a  :=  z  +  1  +  a  -  int (a/4) ; 
end ; 

b  :=  a  +  1524; 

c  :=  i nt ( (b  -  122 . 1) /365 . 25) ; 
d  :=  i nt (365 . 25*c) ; 
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e  :=  trunc((b  -  d)/30.6001); 

Day  :=  b  -  d  -  int(30.6001*e)  +  f; 

if  e  <  13.5  then  Month  :=  e  -  1  else  Month  :=  e  -  13; 
if  Month  >  2.5  then  Year  :=  trunc(c)  -  4716 
else  Year  ;=  trunc(c)  -  4715; 
end; 

{  > 


function  Obliquity  (JD:  real):  real; 

{  This  function  goes  with  the  modified  procedure  Sun,  both  to  be 
used  with  Danby's  GEM  program  to  calculate  the  heliocentric 
ecliptic  orbital  elements  of  a  comet  or  asteroid  referred  to 
the  mean  ecliptic  and  equinox  of  J2000.0. 

Roger  L.  Mansfield,  1998  March  24.  } 

const 

RADIAN  =  57.29577951; 


begin 

Obliquity  :=  23.4392911/RADIAN; 
end ; 

{---- -  - -  — . > 

function  Modf  (X,  Y:  real):  real; 
begi  n 

Modf  :=  Y*Frac(X/Y) 
end ; 

{  } 


Procedure  SUN  (JD:  real;  var  X,  Y,  Z:  real); 

{  Given  a  Julian  Ephemeris  date  JD,  this  procedure  calculates  the 
geocentric  equatorial  cartesian  coordinates  of  the  sun,  in  A.U. 
and  referred  to  the  mean  equator  and  equinox  of  J2000.0. 

Test  data  for  this  procedure  can  be  found  in  the  Astronomical 
Almanac  for  the  Year  1998,  pp.  C20-C23. 

Reference:  Roger  L.  Mansfield,  "Ephemeris  of  a  Comet  via 
Uniform  Path  Mechanics  (UPM),n  a  Mathcad  PLUS  6  worksheet 
at  http://www.mathsoft.com  (Math  in  Action  files.  Astronomy 
and  Navigation  Programs,  September  1997);  see  procedural 
function  HGEO.  } 

const 

RADIAN  =  57.29577951; 

TPI  =  6.283185307; 

SECPDG  =  3600.0; 

OBLQTY  =  23.4392911/RADIAN; 

SECPRV  =  360 . 0*SECPDG ; 

var 

TC,  A,  ECC,  K,  N,  ARGP ,  INC,  LONG,  M,  SIM,  S2M ,  V,  CV,  R, 

XPF ,  YPF ,  CP,  SP,  Cl,  SI,  PX,  PY,  PZ,  QX,  QY,  QZ ,  LM, 

XC,  YC,  ZC,  CO,  SO:  real; 


begi  n 


{  Calculate  Julian  centuries  elapsed  from  input  JD  to  2000 
January  1.5  TT.  } 

TC  :=  (JD  -  2451545. 0)/36525.0; 

{  Calculate  semimajor  axis  and  eccentricity  of  Earth-moon 
barycenter's  orbit  around  sun.  } 

A  :=  1.00000011  -  0 . 00000005*TC ; 

ECC  :=  0.01671022  -  0.00003804*10; 

{  Calculate  mean  motion  of  this  orbit.  Note  that  AA(-3/2)  = 

Exp (-1 . 5*Ln (A) ) .  } 

K  :=  0. 01720209895*Sqrt(l. 00000304) ; 

N  :=  K*Exp(-l. 5*Ln(A) ) ; 

{  Calculate  argument  of  perigee,  inclination,  and  longitude  of 
perihelion.  } 

ARGP  :=  (102.94719  +  1198 . 28*TC/SECPDG) /RADIAN ; 

INC  :=  (0.00005  -  46 . 94*TC/SECPDG) /RADIAN ; 

LONG  :=  (100.46435  +  (1293740.63  +  99 . 0*SECPRV) *TC/SECPDG) /RADIAN ; 

{  Calculate  mean  anomaly.  Calculate  true  anomaly  using  the  equation 
of  the  center  (Moulton,  p.  171).  } 

M  :=  Mod f (LONG -ARGP ,TPI) ; 

SIM  :=  Sin(M); 

S2M  :=  Si n (2 . 0*M) ; 

V  :=  M  +  ECC* (  2 . 0*S1M  +  ECC* (1 . 25*S2M  + 

ECC*(13.0*Sin(3.0*M)  -  3 . 0*S1M) *0 . 083333333  + 

ECC* ( 103 .0*Sin(4.0*M)  -  44 . 0*S2M) *0 . 010416667) ) ; 

{  Calculate  perifocal  cartesian  coordinates  of  Earth.  } 

CV  :=  Cos (V) ; 

R  :  =  A* ( 1 . 0  -  ECC*ECC) / ( 1 . 0  +  ECC*CV) ; 

XPF  :=  R*CV ; 

YPF  :=  R*Si n (V) ; 

{  Calculate  ecliptic  cartesian  coordinates  of  Earth.  } 


CP 

=  Cos (ARGP) 

SP 

=  Si n (ARGP) 

Cl 

=  Cos (INC) ; 

SI 

=  Sin (INC) ; 

PX 

=  CP; 

PY 

=  CI*SP ; 

PZ 

=  SI *SP ; 

QX 

=  -SP; 

QY 

=  Cl *CP ; 

QZ 

=  SI  *CP ; 

X  :=  XPF*PX  +  YPF*QX ; 

Y  :=  XPF*PY  +  YPF*QY ; 

Z  :=  XPF*PZ  +  YPF*QZ ; 

{  Correct  from  Earth-moon  barycenter  to  geocenter  using  the 
mean  orbital  longitude  of  the  moon.  Note  that  this  correction 
is  applied  to  the  ecliptic  cartesian  coordinates.  } 

LM  :=  Modf (218 . 0  +  481268. 0*TC,  360 . 0) /RADIAN ; 

XC  :=  X  -  0 . 0000312*Cos (LM) ; 
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YC  :=  Y  -  0 . 0000312*Si n  (LM) ; 

ZC  :=  Z; 

{  Refer  to  the  mean  equator  and  equinox  of  J2000.0  by  rotating 
about  X-axis  through  angle  OBLQTY,  the  obliquity  of  the 
ecliptic  at  J2000.0.  } 

CO  :=  Cos (OBLQTY); 

SO  :=  Sin(OBLQTY) ; 

X  : =  XC  * 

Y  :=  YC*C0  -  ZC*S0 ; 

Z  :=  YC*S0  +  ZC*C0; 

{  Change  sign  to  convert  heliocentric  equatorial  cartesian  Earth 
vector  to  geocentric  equatorial  cartesian  sun  vector.  } 

X  :=  -X; 

Y  :=  -Y; 

Z  :=  -Z; 

end ; 

{ . - . . - . - . — -  . > 

Procedure  Data; 
begi  n 

Textmode (3) ; 

wr i teln (f i leout ,  '  Input  Data  Echo  Check'); 

for  i  :=  1  to  3  do 
begin 

readln (f i lei n ,  ye,m, fullday, hour , min, sec ,  decdeg , decmi n , dec sec) ; 

d  :=  trunc(fullday) ;  ut  :=  fullday  -  d; 

Date_Jd(ye,  m,  d,  ut,  t i me [ i ] ) ; 
r  a [ i ]  :=  hour  +  (min  +  sec/60)/60; 

r  a [ i ]  :=  r  a  [  i ]  *  15  *  pi / 180 ; 
d e c  [ i ]  :=  decdeg  +  (decmin  +  decsec/60) /60 ; 
dec [ i I.  :=  dec[i]  ^ pi / 180 ; 

{Echo  check  input  data.} 

wri teln  (fi leout ,  ' Observation  ' , i ) ; 

wri teln  (fi leout ,  'Year  =  ' ,  ye,  '  Mon  =  m, 

'  Day  =  1 ,  fullday) ; 

wri teln  (fi leout ,  'HA  =  hour,  'h  ',min,  'm  ',sec,  's'); 

wri teln(f ileout ,  ’DEC  =  decdeg,  'deg  ' .  decmin,  'm  ’, 

decsec ,  ’s'); 

wri teln  (f i leout) ; 

rh[i,l]  :=  cos(ra[i])  *  cos  (dec  [i'] ) ; 
rh  [  i  ,2]  :=  s i n  ( ra  [  i ] )  *  cos(dec[i]); 
rh  [  i  ,  3]  : =  si n  (dec  [  i ] ) ; 

{Direction  cosines.} 

Sun(time[i],  scl[i,l],  scl[i,2],  scl t i . 3] ) ; 
wri teln ; 

end;  {Observation  loop.} 


taul  :=  ( t i me [ 1 ]  -  t i me [ 2 ] )  *  gk; 
tau3  :=  ( t i me [ 3 ]  -  t i me [ 2 ] )  *  gk; 

{(7.3.13).} 
end ; 

{ . .  . — . . . . > 


{Procedure  Data ;■ 
begi  n 

Textmode(3) ; 
for  i  :=  1  to  3  do 
begin 

wri teln( ' Enter  data  for  observation  '  ,i,':'); 
wri teln( '  Time : ' ) ; 

writer  year:  ');  readln(ye); 

writeC  month:  ');  readln(m); 

writer  day:  ');  readln(d).; 

writer  UT:  ');  readln(ut);  writeln; 

Date_Jd(ye,  m,  d,  ut,  t i me [ i ] ) ; 

writelnC  Right  ascension:'); 

writer  hours:  ');  readln  (hour) ; 

writer  minutes:  ');  readln(min); 

writer  seconds:  ');  readln(sec);  writeln; 

r a [ i ]  :=  hour  +  (min  +  sec/60)/60; 

ra[i]  :=  ra [i ] *15*pi/180; 

writelnC  Declination  (degrees  minutes  and  seconds):'); 

writer  degrees:  ’);  readln  (decdeg) ; 

writer  minutes:  ');  readln (degmi n)  ; 

writer  seconds:  ');  readln (degsec) ; 

dec[i]  :=  decdeg  +  (decmin  +  decsec/60) /60 ; 
dec[i]  :=  dec [ i ] *pi /180 ; 
rh  [  i  ,  1]  :=  cos ( ra [ i ] ) *cos (dec  [  i ] ) ; 
rh[i,2]  :=  si n ( ra [ i ] ) *cos (dec [ i ] ) ; 
rh[i  ,3]  :=  sin(dec[i]) ;} 

{Direction  cosines.} 

{Sun (time [i ] ,  scl[i,l],  scl[i,2],  scl[i  ,3] )  ; 
writeln; 

end ;} {Observation  loop.} 

{taul  :=  (time[l]  -  time[2])*gk; 
tau3  :=  (time [3]  -  ti me [2] ) *gk ; } 

{(7.3.13).} 

{end;} 

{ - - - - - - - - 

{Procedure  Data; 
begin 

Textmode(3)  ; } 

{Obs  1} 

{ye  :=  1998;  m  :=  1;  d  :=  21;  ut  :=  0.24164; 

Date_Jd(ye,  m,  d,  ut,  t i me  [  1] ) ; 

hour  :=  3;  min  :=  44;  sec  :=  50.43; 

ra[l]  :=  hour  +  (min  +  sec/60)/60; 

ra [1]  :=  r a [ 1] * 15*pi / 180 ; 

decdeg  :=  42;  decmin  :=  12;  decsec  :=  41.6; 
dec  [  1 ]  :=  decdeg  +  (decmin  +  decsec/60) /60 ; 
dec [ 1]  :=  dec [ 1] *pi /180 ; } 

{Obs  2} 

{ye  :=  1998;  m  :=  2;  d  :=  13;  ut  :=  0.0716; 

Date_Jd(ye,  m,  d,  ut,  time [ 2 ] ) ; 

hour  :=  3;  min  :=  53;  sec  :=  21.41; 

r  a [ 2 ]  :=  hour  +  (min  +  sec/60)/60; 

ra [2]  :=  ra [2] * 15*pi / 180 ; 

decdeg  :=  40;  decmin  :=  32;  decsec  :=  49.2; 
dec[2]  :=  decdeg  +  (decmin  +  decsec/60) /60 ; 


dec  [2]  :=  dec [2] *pi /180 ; } 

{Obs  3} 

{ye  :=  1998;  m  :=  3;  d  :=  13;  ut  :=  0.09222; 

Date_Jd(ye,  m,  d,  ut,  t i me [ 3 ] ) ; 

hour  :=  4;  min  :=  18;  sec  :=  34.54; 

r a [ 3 ]  :=  hour  +  (min  +  sec/60)/60; 
ra [3]  :=  ra [3] *15*pi /180 ; 

decdeg  :=  39;  decmin  : =20 ;  decsec  :=  22.9; 
dec[3]  :=  decdeg  +  (decmin  +  decsec/60) /60; 
dec  [ 3 ]  ;=  dec [3] *pi /180 ; 
for  i  :=  1  to  3  do 
begin 

rh [ i , 1]  :=  cos (ra [ i ] ) *cos (dec [i ] ) ; 
rh[i,2]  :=  si n (ra [ i ] ) *cos (dec [ i ] ) ; 
rh [ i  ,3]  :=  sin(dec[i] ) ; } 

{Direction  cosines.} 

{Sun (time [ i ] ,  scl[i,l]f  scl[i,2],  scl[i,3]); 
wri teln ; 

end;}  {Observation  loop.} 

{taul  :=  (timefl]  -  time [2] ) *gk; 
tau3  :=  ( t i me [ 3 ]  -  ti me [2] ) *gk ; } 

{(7.3.13).} 

{end ; } 

{ - - - - > 

Procedure  Rotation; 

{Resolves  vectors  in  Cunningham's  reference  system.  See  pp  228,  229.} 

var  zl,  z 2 :  real; 

begin 

zl  :=  0; 

for  i  :=  1  to  3  do 
begin 

rm[l , i ]  : =  rh [1 , i ] ; 
zl  :=  zl  +  rh [1 , i ] *rh  [3  ,  i ]  ; 

end;  {Loop  for  xi ,  unit  vector  of  the  first  observation.} 

z2  :=  sqrt (1  -  zl*zl)  ; 

for  i  :  =  1  to  3  do 

rm [ 2 , i ]  :=.(rh[3,i]  -  zl* rh [ 1 , i ] ) /z2 ; 

{Loop  for  eta;  see  (7.3.10).} 

rm [3 , 1]  :=  rm [1 , 2] * rm [2 , 3]  -  rm [ 1 , 3] * rm [ 2 , 2]  ; 
rm  [3 , 2]  :  =  rm  [  1 , 3]  *  rm  [2 , 1]  -  rm  [  1 , 1]  *  rm  [2 , 3]  ; 
rm[3 , 3]  : =  rm[l , 1] * rm [2 , 2]  -  rm[l,2] *rm[2, 1] ; 

{Components  of  zeta;  see  (7.3.11).} 

{Rotate  the  unit  observation  unit  vectors.} 

for  i  :=  1  to  3  do 
begi  n 

rhol[i]  :=  0;  r h o 2 [ i ]  :=  0;  rho3[i]  :=  0; 
end ; 

rhol [ 1]  :=  1; 

{ r hoi [ i ]  are  lambdal,  mul,  nul:  see  (7.3.12).} 
for  i  :=  1  to  3  do 
begi  n 

rho2[l]  :=  rho2[l]  +  rm [ 1 , i ] * rh [2 , i ] ; 
r ho2  [ 2 ]  :=  rho2[2]  +  rm [2 , i ] * rh [2 , i ] ; 
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rho2[3]  :=  rho2[3]  +  rm[3#i]*rh[2.i] ; 

{lambda2,  mu2,  nu2.} 

rho3 [1]  :=  rho3[l]  +  rm [ 1 f i ] *rh [3 t i ] ; 
rho3 [2]  :=  rho3[2]  +  rm [2 . i ] *rh  [3 # i ] ; 

{lambda3  and  mu3;  nu3  =  0.} 
end; 

{Now  rotate  the  solar  coordinates.} 
for  i  :=  1  to  3  do 
begin 

for  j  :=  1  to  3  do 
begin 

sc2[i,j]  :=  0; 
for  k  :=  1  to  3  do 

sc2[i,j]  :=  sc2 [ i ,  j  ]  +  rm[i ,k] *scl[ j ,k]  ; 

end; 

end; 

{Note  that  the  jth  column  of  sc2  represents  the  solar  coordinates 
of  the  jth  time  of  observation.} 
end ; 

{ . ----- . .  . . > 

Function  Q:  real; 

{Calculates  the  function  Q  defined  in  (6.11.20).  Uses  formulas  (6.11.29), 

(6.11.31)  and  (6.11.32).} 

var 

fac,  y,  yl,  qt:  real; 

n:  integer; 

begin 

if  abs(xs)  <  0.1  then 
begi  n 

{Use  hypergeometric  series  (6.11.29).} 
fac  :=  1.2*xs; 
qt  :=  1  +  fac; 

for  n  :=  1  to  10  do 
begi  n 

fac  :=  fac*xs*(3  +  n)/(2.5  +  n)  ; 
qt  :=  qt  +  fac; 

end;  {n  loop  for  hypergeometric  series.} 

Q  :=  qt; 

end 

else 

begin 

if  xs  >  0  then 
begin 

y  :•=  2*xs  -  1; 

yl  :=  Sqrt (xs  -  xs*xs) ; 

Q  :=  (3/16)* (ArcTan(y/Sqrt (1  -  y*y)) 

+  2*y*yl  +  pi /2) / (yl*yl*yl) ; 

{(6.11.31).} 

end 

else 

begi  n 

y  :=  sqrt (xs*xs  -  xs) ; 

Q  :=  (3/16) * (2* ( 1  -  2*xs)*y  -  ln(l  -  2*xs  +  2*y) ) / (y*y*y) ; 

{(6.11.32).} 

end ; 


end ; 
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end ; 


Procedure  Sector_triangle_ratio; 

var  ms,  el,  yg2,  yg3,  den,  dy:  real;  nc:  integer; 


begin 

ms  :=  (gk*dt)*(gk*dt)/ (kay*kay*kay) ; 
el  :=  (rl  +  r2  -  kay) / (2*kay) ; 

{(6.12.11)  and  (6.12.12) .} 

ygi  ;=  l; 

{Initial  guess.} 

{Use  Steffensen's  method.} 
nc  :=  0; 

repeat 

nc  :=  nc  +  1; 

xs  :=  ms/(ygl*ygl)  -  el; 

yg2  :=  1  +  (4*ms*Q)/ (3*ygl*ygl) ; 

{(6.12.14).} 

xs  :=  ms/(yg2*yg2)  -  el; 
yg3  :=  .1  +  (4*ms*Q)/ (3*yg2*yg2) ; 
den  :=  (ygl  -  2*yg2  +  yg3); 
if  abs(den)  >  le-12  then 
begin 

dy  :=  (yg.2  -  ygl)*(yg2  -  ygl)/(den); 
ygl  :=  ygl  -  dy ; 

end; 

until  (abs(dy)  <  le-6)  or  (abs(den)  <  le-12)  or  (nc  =  20); 

if  nc  =  20  then  wri teln (* SECTOR-TRIANGLE  ITERATION  FAILED.*); 
end ; 

{  . - .  . . . . --> 

Function  Atan2(si,  co:  real):  real; 

{Finds  the  angle  between  0  and  2*pi  of  which 

the  sine  is  proportional  to  si  and  the  cosine  is  proportional  to  co.} 
begin 

pi  :=  3.141592653589793; 

if  si  <  0  then 
begin 

if  co  =  0  then  Atan2  :=  1.5*pi 
else 

if  co  >  0  then  Atan2  :=  2 * p i  +  ArcTan (si /co) 
else  Atan2  :=  pi  +  ArcTan (si /co) ; 

end 

else 

begi  n 

if  co  =  0  then  Atan2  :=  pi/2 

else  Atan2  :=  ArcTan (si /co) ; 
if  co  <  0  then  Atan2  :=  pi  +  ArcTan (si /co) ; 
end ; 

end: 

{  . . . . . . . > 

Procedure  Output; 
begin 


wri teln 
wr i teln 
i  nc 

wri teln 
omega 
wri teln 
argp 
wri teln 
wri teln 
wri teln 
Jd_Date 
wri teln 
wri teln 
wri teln 


(AU) 

(unitless) 


=  i nc*180/pi ; 

Inclination  (deg) 
omega*180/pi ; 

Lon  of  A  Node  (deg) 
=  argp*180/pi ; 

Arg  of  Per  (deg) 
n  (deg/day) 

P  (years) 

(tpc,  ye,  m,  d) ; 

('  T,  year 
('  month 

('  day 


.a) ; 

,e) ; 

,  inc) ; 

tomega) ; 

.argp) ; 

.  n) ; 

,  P) ; 

,ye) ; 

.  m) ; 

,d) ; 


wri teln (fi leout , 
wri teln (f i leout , 
wri teln(fileout, 
wri teln(fileout , 
wri teln (f i leout , 
wri teln (f i leout , 
wri teln(f ileout , 
wri teln (fi leout , 
wri teln (fi leout , 
wri teln (fi leout , 
wri teln (fi leout , 
end ; 


Calculated  Orbital- Elements ') ; 


a  (AU) 

=  1 

.a); 

e  (unitless) 

=  ' 

.  e) ; 

Inclination  (deg) 

=  1 

,  inc) ; 

Long  of  A  Node  (deg) 

=  ' 

.omega) 

Arg  of  Per  (deg) 

=  1 

•argp) ; 

n  (deg/day) 

=  ' 

,  n) ; 

P  (years) 

=  1  , 

.  P): 

T,  year 

=  '  , 

,  ye) ; 

month 

=  * 

,  m) ; 

day 

=  1 

,  d); 

Procedure  Coordi nates_elements ; 
var 

mu,  hx,  hy,  hz,  h,  ex,  ey,  ez,  u,  s,  c,  ean:  real; 
begin 

mu  :=  gk*gk; 

hx  :=  y*zv  -  z*yv; 

hy  :=  z*xv  -  x*zv; 

hz  :=  x*yv  -  y*xv; 

h  :=  Sqrt(hx*hx  +  hy*hy  +  hz*hz)  ; 
r  :=  Sqrt(x*x  +  y*y  +  z*z); 
ex  :=  (yv*hz  -  zv*hy)/mu  -  x/r; 
ey  :=• (zv*hx  -  xv*hz)/mu  -  y/r; 
ez  :=  (xv*hy  -  yv*hx)/mu  -  z/r; 
u  :=  x*xv  +  y*yv  +  z*zv; 

s  :=  hx;  c  :=  -  hy; 

omega  :=  Atan2(s,c) ; 

s  :=  Sqrt(hx*hx  +  hy*hy) ;  c  :=  hz; 

inc  :=  Atan2(s ,c) ; 

s  ;=  ez;  c  :=  (-  ex*hy  +  ey*hx)/h; 

argp  :=  Atan2(s ,c) ; 

e  :=  Sqrt(ex*ex  +  ey*ey  +  ez*ez) ; 
a  :=  l/(2/r  -  (xv*xv  +  yv*yv  +  zv*zv)/mu); 
n  :=  gk* (1/Sqrt (a*a*a) ) ; 
n  :=  n  *  180.0/pi ; 
p  :=  360/ n*  (1/36-5.25)  ; 

if  a  >  0  then 
begin 

s  :=  u/Sqrt (a*mu) ;  c  :=  1  -  r/a; 
ean  : =  Atan2 (s , c) ; 


tpc  :=  t  -  (ean  -  s) *Sqrt (a*a*a/mu) ; 

end 

else 

begin 

z  (1  -  r/a)/e; 

ean  :=  ln(z  +  Sqrt(z*z  -  1)); 

if  u  <  0  then  ean  :=  -  ean; 

tpc:=  t  -  (e*(exp(ean)  -  exp(-  ean))/2  -  ean)*Sqrt(-  a*a*a)/gk; 
end ; 

Output ; 
end; 

{ - - - > 

begin  {Main  program.} 

assign(f ilein,  ’input. txt') ;  reset (fi lei n) ; 

assign(fileout,  ' output . txt * ) ;  rewr i te (f i leout) ; 

Data; 

Rotation ; 

{Ini ti ali ze  . } 

cl  :=  tau3/(tau3  -  taul) ; 
c3  :=  -  taul/(tau3  -  taul); 

{(7.3.14).} 

for  i  ;=  1  to  3  do  gd[i]  :=  0; 
writeln; 

COUNT  :=  0; 

repeat 

COUNT  :=  COUNT  +  1; 
wri teln(  'cl  =  '  ,cl)  ; 
wri teln(  'c3  =  '  ,c3)  ; 

temp  :=  (-  cl*sc2[3,l]  +  sc2[3,2]  -  c3*sc2  [3 f 3] ) /rho2  [3] ; 
change  :=  abs(temp  -  gd [ 2 ] ) ; 
gd [2]  :=  abs(temp) ; 

{(7.3.15.} 

temp  :=  (gd [2] * rho2 [2]  +  cl*sc2[2,l]  -  sc2[2,2]  +  c3*sc2[2,3]) 

/ (c3*  rho3 [2] ) ; 

change  :=  change  +  abs(temp  -  gd [3] ) ; 
gd [ 3 ]  :=  abs(temp) ; 

{(7.3.16).} 

temp  :=  (gd [2] *rho2 [1]  -  c3*gd [3] *rho3 [1] 

+  cl*sc2 [1,1]  -  sc2  [  1 , 2]  +  c3*sc2[l,3])/(cl*rhol[l]); 
change  :=  change  +  abs(temp  -  gd [ 1 ] ) ; 
gd [1]  :=  abs(temp) ; 

{(7.3.17)  .} 

for  i  :=  1  to  3  do  wri teln ( ' rho ' , i , '  =  '.gd[i]); 

if  (change  >  0.001)  then 
begi  n 

{Find  heliocentric  coordinates.} 
for  i  :=  1  to  3  do 
begi  n 

he [i , 1]  :=  gd [ 1] * rhol [ i ]  -  sc2 [ i , 1] ; 
he [i , 2]  :=  gd [2] * rho2 [i ]  -  sc2 [ i , 2]  ; 
he [ i ,3]  :=  gd [3] *rho3  [i ]  -  sc2 [ i , 3]  ; 
end ; 


{The  jth  column  of  he  represents  the  jth  observation.} 

{Now  find  the  sector-tri angle  ratios.} 

{Find  Yl} 

rl  :=  0;  r2  :=  0;  kay  :=  0; 

for  i  :=  1  to  3  do 
begin 

rl  :=  rl  +  he [i , 2] *hc [i  , 2]  ; 
r2  :=  r2  +  he [ i ,3] *hc[i ,3]  ; 
kay  :  =  kay  +  he [i , 2] *hc [i  , 3] ; 
end ; 

rl  :=  sqrt (rl) ; 

r2  :=  sqrt ( r 2 )  ; 

kay  :=  sqrt(2*kay  +  2*rl*r2); 

dt  :=  time [3]  -  time [2] ; 

Sector_tri angle_ratio; 
yl  :=  ygi; 
writelnCYl  =  '  .  yl) ; 

{Find  Y2 } 

rl  :=  0;  r2  :=.0;  kay  :=  0; 

for  i  :=  1  to  3  do 
begin 

rl  :=  rl  +  he [ i ,1] *hc[i ,1] ; 
r2  : =  r2  +  hc[i .3] *hc[i ,3] ; 
kay  :=  kay  +  he [i , 1] *hc  [i ,  3]  ; 
end ; 

rl  :=  sqrt  (rl) ; 

r2  :=  sqrt (r2)  ; 

kay  :=  sqrt(2*kay  +  2*rl*r2); 

dt  :=  time[3]  -  time[l] ; 

5ector_tri angle_ratio; 

y2  :=  ygl; 

wr i teln ( 1 Y2  =  \  y2) ; 


{Find  Y3} 

rl  :=  0;  r2  :=  0;  kay  :=  0; 

for  i  :=  1  to  3  do 
begi  n 

rl  :=  rl  +  he [i ,l]*hc[i ,1]  ; 

r2  :=  r2  +  he [ i ,2]*hc[i  ,2]  ; 

kay  :=  kay  +  he [ i  , 1] *hc [ i , 2]  ; 

end ; 

rl  :=  sqrt(rl) ; 

r 2  :=  sqrt (r2)  ; 

kay  :=  sqrt(2*kay  +  2*rl*r2); 

dt  :=  time[2]  -  time[l] ; 

Seetor_t ri angle_ratio ; 

y3  :=  ygi; 

wri teln ( ’ Y3  =  ’ ,y3) ; 

{Update  guess  for  cl  and  c3} 

cl  :=  (y2/yl) * ( tau3) / (tau3  -  taul); 
c3  :=  (y2/y3)*(  -  taul)/(tau3  -  taul); 

end;  {Loop  for  iterating  on  cl  and  c3.} 

wri teln (' change  =  1 , change); 


writeln; 

wr  i teln ( ' COUNT  =  ’ ,C0UNT:3) ; 

write  (‘CONTINUE  (1)  OR  STOP  (0)  ?  “); 

readln  (CONTINUE); 

writeln; 

until  (change  <  0.001)  or  (keypressed)  or  (COUNT  =  100)  or  (CONTINUE 
{The  'small  quantity'  can  be  changed  at  your  discretion.} 

{Find  heliocentric  equatorial  coordinates  for  times  1  and  2.} 

for  i  :=  1  to  3  do 
begin 

x 1 [ i ]  :=  gd [1] *rh [1 , i ]  -  scl[l,i]  ; 
x2[i]  :=  gd [2] *rh  [2 ,  i ]  -  scl [2 , i ] ; 
end ; 

r  sqrt (xl  [1] *xl  [1]  +  xl[2]*xl[2]  +  xl [3] *xl [3] ) ; 
f  :=  1  -  2*gk*gk*dt*dt/(kay*kay*ygl*ygl*r) ; 
g  :=  dt/ygl ; 

{Find  velocity  components  at  time  1.} 

for  i  :=  1  to  3  do 

x 2 [ i ]  :=  (x2 [ i ]  -  f*xl[i])/g; 

{Rotate  to  ecliptic  coordinates.} 
obi  :=  Obliqui ty (time [2] )  ; 
x  :  =  xl [1] ; 

y  :=  xl [2] *cos (obi)  +  xl [3] *si n (obi) ; 
z  :=  -  xl [2] *si n (obi)  +  xl [3] *cos (obi ) ; 
xv  :=  x2  [  1] ; 

yv  :=  x  2 [ 2 ] *cos (obi)  +  x2 [3] *sin(obl) ; 
zv  :=  -  x2  [2] *sin(obl)  +  x2 [3] *cos (obi) ; 

{Find  elements.} 
t  :=  time [1] ; 

Coordi nates_elements ; 

close(filein) ; 
close (fi leout)  ; 
end . 


APPENDIX  E 


GEM  Output  Data 


TEST  CASE  1 
Input  Data  Echo  Check 
Observation  1 

Year  =  1998  Mon  =  1  Day  =  2.1241640000E+01 

HA  =  3.0000000000E+00h  4.4000000000E+01m  5.0430000000E+01S 

DEC  =  4.2000000000E+01  deg  1.2000000000E+01m  4.1600000000E+01s 

Observation  2 

Year  =  1998  Mon  =  2  Day  =  1.306861 0000E+01 

HA  =  3.0000000000E+00h  5.3000000000E+01m  2.14100Q0000E+01S 

DEC  =  4.0000000000E+01  deg  3.2000000000E+01m  4.9200000000E+01s 


Observation  3 

Year  =  1998  Mon  =  3  Day  =  1 .3092220000E+01 

HA  =  4.0000000000E+00h  1.8000000000E+01m  3.4540000000E+01s 

DEC  =  3.9000000000E+01  deg  2.0000000000E+01m  2.2900000000E+01S 


Calculated  Orbital  Elements 


a  (AU) 
e  (unitless) 

Inclination  (deg) 

Long  of  A  Node  (deg) 
Arg  of  Per  (deg) 
n  (deg/day) 

P  (years) 

T, year 
month 
day 


=  3.21 82660546E+00 
=  1.5972662894E-01 
=  1. 806201 2065E+01 
=  1.4028081 831 E+00 
=  3.2927072274E+02 
=  1.7071479356E-01 
=  5.7735259071  E+00 
=  1996 
=  8 

=  6.872951 5076E+00 
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TEST  CASE  2 
Input  Data  Echo  Check 
Observation  1 

Year  =  1998  Mon  =  1  Day  =  2.2217710000E+01 

HA  =  3.0000000000E+00h  4.4000000000E+01m  5.5410000000E+01s 

DEC  =  4.2000000000E+01  deg  7.0000000000E+00m  4.2600000000E+01s 

Observation  2 

Year  =  1998  Mon  =  2  Day  =  1.3071 600000E+01 

HA  =  3.0000000000E+00h  5.3000000000E+01  m  2.1530000000E+01S 

DEC  =  4.0000000000E+01  deg  3.2OOO0OOOOOE+O1m  4.8700000000E+01s 


Observation  3 

Year  =  1998  Mon  =  3  Day  =  1.3094720000E+01 

HA  =  4.0000000000E+00h  1.8000000000E+01m  3.471 0000000E+0 Is 

DEC  =  3.9000000000E+01  deg  2.0000000000E+01m  2.2400000000E+01S 


Calculated  Orbital  Elements 


a  (AU) 
e  (unitless) 

Inclination  (deg) 

Long  of  A  Node  (deg) 
Arg  of  Per  (deg) 
n  (deg/day) 

P  (years) 

T, year 
month 
day 


=  3.2101537365E+00 
=  1. 63291 92875E-01 
=  1. 80631 42247E+01 
=  1. 483378701 7E+00 
=  3.2862022751 E+02 
=  1.71 36231 735E-01 
=  5.751 7095859E+00 
=  1996 
=  8 

=  7.6320877075E+00 
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TEST  CASE  3 


Input  Data  Echo  Check 
Observation  1 

Year  =  1998  Mon  =  1  Day  =  2.721 5900000 E+01 

HA  =  3.0000000000E+00h  4.5000000000E+01m  4.6100000000E+01S 

DEC  =  4. 1 000000000E+01  deg  4.3000000000E+01m  2.5000000000E+00s 


Observation  2 

Year  =  1998  Mon  =  2  Day=  1 .3090280000E+01 

HA  =  3.0000000000E+00h  5.3000000000E+01m  2.2220000000E+01s 

DEC  =  4.0000000000E+01  deg  3.2000000000E+01m  4.4800000000E+01s 


Observation  3 

Year  =  1998  Mon  =  3  Day  =  1 .3096840000E+01 

HA  =  4.0000000000E+00h  1 ,8000000000E+01m  3.4840000000E+01S 

DEC  =  3.9000000000E+01  deg  2.0000000000E+01m  2.2200000000E+01s 


Calculated  Orbital  Elements 


a  (All) 
e  (unitless) 

Inclination  (deg) 

Long  of  A  Node  (deg) 
Arg  of  Per  (deg) 
n  (deg/day) 

P  (years) 

T, year 
month 
day 


=  3.2268883584E+00 
=  1.5592254265E-01 
=  1 .8061 1 30253E+01 
=  1. 36481 44047E+00 
=  3.29910561 36E+02 
=  1.70031 021 33E-01 
=  5.7967438862E+00 
=  1996 
=  8 

=  5.9996833801  E+00 
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TEST  CASE  4 


Input  Data  Echo  Check 
Observation  1 

Year  =  1998  Mon  =  1  Day  =  2.7281 770000E+01 

HA  =  3.0000000000E+00h  4.5000000000E+01m  4.7030000000E+01S 

DEC  =  4.1000000000E+01deg  4.2000000000E+01m  4.3500000000E+01s 


Observation  2 

Year  =  1998  Mon  =  2  Day  =  1 .3092640000E+01 

HA  =  3.0000000000E+00h  5.3000000000E+01m  2.2330000000E+01S 

DEC  =  4.0000000000E+01  deg  3.2000000000E+01m  4.4000000000E+01S 

Observation  3 

Year  =  1998  Mon  =  3  Day=  1.3132670000E+01 

HA  =  4.0000000000E+00h  1.8000000000E+01m  3.7230000000E+01S 

DEC  =  3.9000000000E+01  deg  2.0000000000E+01m  1.8700000000E+01s 


Calculated  Orbital  Elements 


a(AU) 
e  (unitless) 

Inclination  (deg) 

Long  of  A  Node  (deg) 
Arg  of  Per  (deg) 
n  (deg/day) 

P  (years) 

T,  year 
month 
day 


=  3.2727041 601 E+00 
=  1.2374729827E-01 
=  1.8049346092E+01 
=  9.5408620780E-01 
=  3.3445223269E+02 
=  1.6647305562E-01 
=  5.9206355028E+00 
=  1996 
=  7 

=  2.9782054901 E+01 


TEST  CASE  5 


Input  Data  Echo  Check 
Observation  1 

Year  =  1998  Mon  =  1  Day  =  2.821 3840000E+01 

HA  =  3.0000000000E+00h  4.6000000000E+01m  1.1700000000E+00s 

DEC  =  4.1000000000E+01deg  3.8000000000E+01m  1.8700000000E+01s 

Observation  2 

Year  =  1998  Mon  =  2  Day  =  1.3123760000E+01 

HA  =  3.0000000000E+C>0h  5.3000000000E+01 m  2.3490000000E+01S 

DEC  =  4.0000000000E+01  deg  3.2000000000E+01m  3.7600000000E+01s 


Observation  3 

Year  =  1998  Mon  =  3  Day  =  1.31 38630000E+01 

HA  =  4.00000000()0E+00h  1 .8000000000E+01  m  3.7630000000E+01s 

DEC  =  3.9000000000E+01  deg  2.0000000000E+01m  1.7600000000E+01s 


Calculated  Orbital  Elements 


a  (AU) 
e  (unitless) 

Inclination  (deg) 

Long  of  A  Node  (deg) 
Arg  of  Per  (deg) 
n  (deg/day) 

P  (years) 

T,  year 
month 
day 


3.20711 821 32E+00 
1. 6607.1 23946E-01 
1.8061309603E+01 
1 .5977 1 244 1 4E+00 
3.2851 102974E+02 
1.7160566543E-01 
5.7435532848E+00 
1996 
8 

1 .01 76685333E+01 
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TEST  CASE  6 
Input  Data  Echo  Check 
Observation  1 

Year  =  1998  Mon  =  1  Day  =  2.8270690000E+01 

HA  =  3.0000000000E+00h  4.6000000000E+01m  2.0200000000E+O0s 

DEC  =  4.1 000000000E+01  deg  3.8000000000E+01m  2.5000000000E+O0s 

Observation  2 

Year  =  1998  Mon  =  2  Day  =  1.3126380000E+01 

HA  =  3.0000000000E+C)0h  5.3000000000E+01m  2.3590000000E+01S 

DEC  =  4.0000000000E+01  deg  3.2000000000E+01m  3.7100000000E+01s 

Observation  3 

Year  =  1998  Mon  =  3  Day  =  1.3140650000E+01 

HA  =  4.0000000000E+00h  1 .8000000000E+01  m  3.7760000000E+01s 

DEC  =  3.9000000000E+01  deg  2.0000000000E+01m  1.7300000000E+01S 

Calculated  Orbital  Elements 
a  (AU) 
e  (unitless) 

Inclination  (deg) 

Long  of  A  Node  (deg) 

Arg  of  Per  (deg) 
n  (deg/day) 

P  (years) 

T,  year 
month 
day 


=  3.2503598737E+00 
=  1. 431401 5144E-01 
=  1. 80535901 08E+01 
=  1 .1 920248958E+00 
=  3.321 8643621 E+02 
=  1.6819260693E-01 
=  5.8601046821  E+00 
=  1996 
=  8 

=  4.79381 17981  E+00 


APPENDIX  F 


Modified  ORBDET  Program 


(* - - - - 

(*  ORBDET 

(*  Gaussian  orbit  determination  from  three  observations 

(*  using  the  abbreviated  method  of  Bucerius 

(*  version  93/07/01 

(*-- . . - . . . . 


PROGRAM  ORBDET (INPUT , OUTPUT, ORBINP , 0RB0UT) ; 


USES  MATLIB,  PNULIB,  SPHLIB,  SUNLIB,  TIMLIB,  KEPLIB; 


TYPE  CHAR80  =  ARRAY[1..80]  OF  CHAR; 


VAR  TEQX 

TP , Q , ECC , INC , LAN , AOP 
J  DO 

RSUN ,  E 
HEADER 

ORBINP , ORBOUT 
MA 


:  REAL; 

REAL; 

:  REAL3 ; 

:  MAT3X; 

:  CHAR80 ; 

TEXT; 

:  REAL; 


*) 

*) 

*) 

*) 

*) 

*) 


(* - - - - - - - - - - *) 

(*  START:  reads  the  input  file  and  preprocesses  the  observational  data  *) 

(*  *) 

(*  output:  *) 

(*  RSUN:  matrix  of  three  Sun  position  vectors  in  ecliptic  coordinates  *) 

(*  E:  matrix  of  three  observation  direction  unit  vectors  *) 

(*  JD:  julian  date  of  the  three  observation  times  *) 

(*  TEQX:  equinox  of  RSUN  and  E  (in  Julian  centuries  since  J2000)  *) 

(* - - - - *) 


PROCEDURE  START  (VAR  HEADER:  CHAR80 ; 

VAR  RSUN , E :  MAT3X ;  VAR  JDO:  REAL3;  VAR  TEQX:  REAL); 


VAR  DAY „ MONTH , YEAR , D , M , I  :  INTEGER; 
UT,S, DUMMY 
EQXO , EQX , TEQXO 

ls,bs,rs,lp,bp,ra,dec,t 

A,  AS 

ORBINP 


REAL; 

REAL; 

REAL3 ; 

:  REAL33 ; 
:  TEXT; 


BEGIN 

(*  open  input  file  *) 

(*  RESET (ORBINP) ;  *)  (*  Standard  Pascal  *) 

ASSIGN (ORBINP , 'ORBINP . DAT ' ) ;  RESET (ORBINP) ;  (*  Turbo  Pascal  *) 


(*  read  data  from  file  ORBINP  *) 

FOR  I :=1  TO  80  DO  (*  header  *) 

IF  NOT (EOLN (ORBINP) )  THEN  READ (ORBINP , HEADER [ I ] )  ELSE  HEADER [ I ] : = 1 
READLN (ORBINP) ; 

FOR  I  :=  1  TO  3  DO 


(*  3  observations  *) 
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BEGIN 

READ  (ORBINP, YEAR, MONTH, DAY, UT) ;  (*  date  *) 

READ  (ORBINP, D, M, S) ;  DDD(D , M , S , RA [ I] ) ;  (*RA  *) 

READLN(ORBINP,D,M,S) ;  DDD(D, M , S , DEC [I] ) ;  (*  Dec  *) 

RA [ I] : =RA [ I] *15 . 0 ; 

J DO [ I ]  :=  2400000. 5+MJD( DAY, MONTH, YEAR, UT) ; 

T[I]  :=  ( J  DO [I] -2451545 .0) /36525 .0; 

END; 

WRITELN ; 

READLN(ORBINP.EQXO) ;  TEQXO : = ( EQXO- 2000 . 0) / 100 . 0 ;  (*  equinox  *) 

(*  desired  equinox  of  the  orbital  elements  *) 

READ(ORBINP , EQX  );  TEQX  :=(EQX  -2000.0)/100.0; 

(*  calculate  initial  data  of  the  orbit  determination  *) 

PMATECL (TEQXO, TEQX, A) ; 

FOR  I  :=  1  TO  3  DO 
BEGIN 

CART  (1.0, DEC [I] ,RA[I] ,E[I,X] , E [ I . Y] , E [I . Z] ) ; 

EQUECL  (  TEQXO  , E [I , X] , E [ I , Y] , E [I , Z] ) ; 

PRECART (  A  , E [I , X] , E [ I , Y] , E [I , Z] )  ; 

POLAR  (E [I ,X] ,E[I,Y] ,E[I,Z] , DUMMY, BP [I] , LP [I] ) ; 

PMATECL (  T [ I] .TEQX, AS) ; 

SUN200  (T[ I ] , LS [I] , BS [I ]  ,  RS [ I ]  )  ; 

CART  (RS [I] , BS [I] , LS [ I ] ,RSUN[I,X] .RSUNfl.Y] , RSUN [ I , Z] ) ; 

PRECART (  AS  , RSUN [I , X] , RSUN [I , Y] , RSUN [ I , Z] ) ; 

END; 

(*  open  file  for  writing  *) 

(*  REWRITE (ORBOUT) ;  *)  (*  Standard  Pascal  *) 

ASSIGN (ORBOUT, 1 ORBOUT.DAT' ) ;  REWRITE (ORBOUT) ;  (*  Turbo  Pascal  *) 

WRITELNC  ORBDET :  orbit  determination  from  three  observations  ’); 

WRITELNC  version  93/07/01  '); 

WRITELN ( '  (c)  1993  Thomas  Pfleger,  Oliver  Montenbruck  ’); 

WRITELN;  WRITELN; 

WRITELNC  Summary  of  orbit  determination  '); 

WRITELN; 

WRITE  ('  ’);  FOR  I:=l  TO  78  DO  WRITE (HEADER [I] ) ;  WRITELN; 

WRITELN; 

WRITELNC  Initial  data  (ecliptic  geocentric  coordinates  (in  deg))'); 
WRITELN; 

WRITELNC  Julian  Date  ' ,  JDO [1] : 12 ; 2 . JD0[2] : 12 : 2 , JDO [3] : 12 : 2) ; 

WRITELNC  Solar  longitude  ',  LS [ 1] : 12 : 2 ,  LS[2]:12:2,  LS [3]  :  12 : 2) ; 
WRITELNC  Planet/Comet  Longitude' . LP [ 1 ] :9:2,  LP [2] : 12 : 2 ,  LP[3]:12:2); 
WRITELNC  Planet/Comet  Latitude  ',BP[1]:9:2,  BP [ 2 ] : 12 : 2 .  BP [3] : 12 : 2) ; 
WRITELN;  WRITELN; 

WRITELN(0RB0UT, '  ORBDET:  orbit  determination  from  three  observations  '); 
WRITELN(0RB0UT , '  version  93/07/01  '); 

WRITELN(0RB0UT, '  (c)  1993  Thomas  Pfleger,  Oliver  Montenbruck  '); 

WRITELN(ORBOUT) ;  WRITELN (ORBOUT) ; 

WRITELN(0RB0UT, '  Summary  of  orbit  determination  '); 

WRITELN (ORBOUT) ; 

WRITE  (ORBOUT,'  ');  FOR  I:=l  TO  78  DO  WRITE(HEADER[I] ) ;  WRITELN; 
WRITELN(ORBOUT) ; 

WRITELN(0RB0UT, '  Initial  data  (ecliptic  geocentric  coordinates  (in  deg))'); 
WRITELN(ORBOUT) ; 

WRITELN (ORBOUT , '  Julian  Date 

J  DO [ 1] : 12: 2, JDO [2] : 12: 2, JDO [3] ; 12 : 2) ; 
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WRITELN(ORBOUT, '  Solar  longitude  LS [1] : 12 : 2 ,  LS [2] : 12 : 2 , 

LS [3] : 12 : 2) ; 

WRITELN(0RB0UT, '  Planet/Comet  Longitude' ,LP[1]  :9:2,  LP [ 2] : 12 : 2 , 

LP [3] : 12 : 2) ; 

WRITELN(ORBOUT, 1  Planet/Comet  Latitude  ' , BP [ 1] : 9 : 2 .  BP [2] : 12 : 2 . 

BP [ 3 ] : 12 : 2) ; 

WRITELN(ORBOUT) ;  WRITELN (ORBOUT) ; 


END; 


(*-  *> 
(*  DUMPELEM:  output  of  orbital  elements  (screen)  *) 

(* - - - - - - - - - - . . *) 


PROCEDURE  DUMPELEM(TP, Q, ECC , INC , LAN .AOP ,TEQX: REAL) ; 
CONST  KGAUSS  =  0.01720209895;  PI  =  3.141592654; 
VAR  DAY, MONTH. YEAR  :  INTEGER; 

MODJD.UT  :  REAL; 

AX, MM, MA, PER  :  REAL; 

BEGIN 


MODJD  :  =  TP*36525 .0  +  51544.5; 

CALDAT(  MODJD,  DAY , MONTH , YEAR , UT) ; 

AX  :=  Q/ (1.0- ECC) ; 

MM  :=  KGAUSS  /  SQRT (ABS (AX*AX*AX) ) ; 

MA  ;=  MM  *  (2450800 . 5- (M0DJD+2400000 . 5) ) ; 

MM  :=  MM  *  180.0/PI; 

MA  :=  MA  *  180.0/PI; 

PER  :=  360.0/MM  *  (1/365.25) ; 

WRITELNC  Orbital  elements', 

'  (Equinox  ' , ' J ' , 100 . 0*TEQX+2000 . 0 : 8 : 2 , ' ) ' ) ; 

WRITELN; 


WRITELNC 

Perihelion  date 

tp 

• 

YEAR: 4, 1 / 1 .MONTH: 2, . 

DAY:  2, 

UT: 8 : 4 , 'h1 . 

i 

(JD1 , MODJ  D+2400000 . 5 

:  11 : 2 , 

')'): 

WRITELNC 

Perihelion  distance 

q  [AU] 

\  Q: 12 : 6) ; 

WRITELNC 

Semi-major  axis 

a  [AU] 

1 .  Q/C1-ECC) 

12:6); 

WRITELNC 

Eccentricity 

e 

V  ECC : 12 : 6) 

WRITELNC 

Inclination 

i 

\  INC : 10 : 4 , 

degrees  1 ) ; 

WRITELNC 

Ascending  node 

Omega 

1 ,  LAN: 10:4 C 

degrees  1 ) ; 

WRITELNC1 

Long,  of  perihelion 

Pi 

A0P+LAN:10:4, 1  degrees1) 

WRITELNC1 

Arg.  of  perihelion 

omega 

\  AOP : 10 : 4 , 

degrees  1 ) ; 

WRITELNC 

Mean  motion 

n 

\  MM : 12 : 6 , 1 

degrees/day 1 ) ; 

WRITELNC 

Mean  anomaly 

M 

\  MA: 12 : 6 , 1 

degrees  1 ) ; 

WRITELNC 

WRITELN; 

END; 

Orbital  period 

P 

\  PER: 10:4, 

years  1 ) ; 

(*--- . - . - .  - .  - - - *> 

(*  SAVEELEM:  output  of  orbital  elements  (file)  *) 

(*-- - - - - - - *) 


PROCEDURE  SAVEELEM (TP , Q ,  ECC , INC , LAN , AOP .TEQX: REAL ; 
HEADER:  CHAR80) ; 

CONST  KGAUSS  =  0.01720209895;  PI  =  3.141592654; 
VAR  I, DAY. MONTH, YEAR  :  INTEGER; 

MODJD.UT  :  REAL; 

AX, MM, MA, PER  :  REAL; 

BEGIN 

MODJD  :=  TP*36525 . 0  +  51544.5; 

CALDAT (  MODJD,  DAY , MONTH , YEAR , UT) ; 
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AX  :=  Q/ (1 .O-ECC) ; 

MM  :=  KGAUSS  /  SQRT (ABS (AX*AX*AX) ) ; 

MA  :=  MM  *  (2450800.5- (M0DJD+2400000 . 5)) ; 

MM  :=  MM  *  180.0/PI; 

MA  :=  MA  *  180.0/PI; 

PER  :=  360.0/MM  *  (1/365.25)  ; 

WRITE  (0RB0UT, YEAR: 5 .MONTH : 3 , (DAY+UT/24 .0) ;7: 3 , ' !  '  :6); 
WRITELN(ORBOUT, 1  perihelion  time  TO  (y  m  d.d)  =  JD  ' , 
(M0DJD+2400000 . 5) ; 12 : 3) ; 

WRITELN(ORBOUT.  Q  : 12 : 6 .  ’  !  ’  :  9,'  q  (a  =' ,Q/(1-ECC) :10:6, ’  )'); 
WRITELN(ORBOUT, ECC: 12:6, ' ! ' :  9,  '  e  '); 

WRITE LN (ORBOUT, INC : 10: 4 , ’ ! ' : 11 , '  i  '); 

WRITE LN (ORBOUT, LAN : 10:4 , ' ! 1 : 11 , '  long. asc . node  '); 

WRITE LN (ORBOUT ,AOP : 10:4 , 1 ! 1 : 11 . 

1  arg.perih.  (  long. per.  =  1 ,A0P+LAN:9:4, 1  )’); 
WRITELN(0RB0UT,MM:12:6. * ! ‘ :  9,'  n  '); 

WRITELN(0RB0UT, MA:  12 : 6 ,  '  !.'  :  9,'  M  '); 

WRITE LN (ORBOUT, PER: 10:4 , ' ! ' : 11 , '  P  '); 

WRITELN(ORBOUT, TEQX*100. 0+2000. 0:8:2, '!' :13, ’  equinox  ( J )  ’ ) ; 
WRITE  (ORBOUT,'!  '); 

FOR  I : =1  TO  78  DO  WRITE (ORBOUT, HEADER [ I] ) ; 

RESET (ORBOUT) ;  (*  close  file  *) 


END; 


(*  RETARD:  light-time  correction  *) 
(*  JDO:  times  of  observation  (tl ' , t2 ' , t3 ' )  (Julian  Date)  *) 
(*  RHO:  three  geocentric  distances  (in  AU)  .  *) 
(*  JD:  times  of  light  emittance  (tl,t2,t3)  (Julian  Date)  *) 
(*  TAU:  scaled  time  differences  *) 


(* - -  - - - - - - - 

PROCEDURE  RETARD  (  JDO, RHO:  REAL3 ;  VAR  JD.TAU:  REAL3) ; 

CONST  KGAUSS  =  0.01720209895;  A  =  0.00578; 

VAR  I:  INTEGER; 

BEGIN 

FOR  I : =1  TO  3  DO  JD[I] : = J DO [ I ] - A* RHO [ I ] ; 

TAU [ 1]  :=  KGAU  SS*(JD[3] -JD[2] ) ;  TAU[2]  :=  KGAUSS* (JD[3]-JD[1]); 
TAU [3 ]  :=  KGAU  SS*(JD[2] - J  D [ 1 ] ) ; 

END; 


(* . . 

(*  GAUSS: 

(* 

(*  RSUN : 
(*  E  : 
(*  JDO  : 
(.*  TP  : 

(*  Q  : 
(*  ECC  : 
(*  INC  : 
(*  LAN  : 
(*  AOP  : 
(* - 


iteration  of  the  abbreviated  Gauss  method 

three  vectors  of  geocentric  Sun  positions 

three  unit  vectors  of  geocentric  observation  directions 

three  observation  times  (Julian  Date) 

time  of  perihelion  passage  (Julian  centuries  since  J2000) 

perihelion  distance 

eccentricity 

inclination 

longitude  of  the  ascending  node 
argument  of  perihelion 


PROCEDURE  GAUSS  (  RSUN.E:  MAT3X ;  JDO : REAL3 ; 

VAR  TP. Q, ECC, INC, LAN, AOP;  REAL  ); 

CONST  EPS_RH0  =1 . 0E- 8 ; 


*) 

*) 

*) 

•) 

*) 

•) 

*) 

*) 

*) 

*) 

*) 

*) 

*) 


VAR  I.J 


:  INTEGER; 


s  :  INDEX; 

RHOOLD ,  DET  :  REAL; 

JD ,  RHO , N , TAU , ETA  :  REAL3 ; 

DI  :  VECTOR; 

RPL  :  MAT3X; 

DD  :  REAL33 ; 

BEGIN 

(*  calculate  initial  approximations  of  nl  and  n3  *) 

N [ 1]  :=  (JDO [3] - JDO [2] )  /  (JDO[3J - JDO [1] ) ;  N[2]  :=  -1.0; 

N [3]  :=  (JDO [2] -JDO [1] )  /  (JD0[3] -JD0[1] ) ; 

(*  calculate  matrix  D  and  its  determinant  (det(D)  =  e3.d3)  *) 

CROSS (E [2] . E [3] ,DI) ;  FOR  J : =1  TO  3  DO  DD [ 1 , J ] : =DOT(DI , RSUN [ J] ) ; 

CROSS (E [3] , E [1] ,  DI) ;  FOR  J : =1  TO  3  DO  DD [2 , J ] : =DOT(DI , RSUN [ J ] ) ; 

CROSS (E [ 1] , E [2] ,  DI )  ;  FOR  J : =1  TO  3  DO  DD [3 , J ] :=DOT(DI , RSUN [ J] ) ; 

DET  :=  DOT (E [3] ,DI) ; 

WRITELN ;  WRITELN ( 1  Iteration  of  the  geocentric  distances  rho  [AU]  ') 
WRITELN ; 

RHO [2]  :=  0; 

(*  Iterate  until  distance  rho[2]  does  not  change  any  more  *) 

RHO  [2]  :=  0; 

REPEAT 

RHOOLD  :=  RHO [2] ; 

(*  geocentric  distance  rho  from  nl  and  n3  *) 

FOR  I  :=  1  TO  3  DO 

RHO [ I ] :  =  (  N [ 1] *  DD [ 1 , 1 ]  -  DD [ 1 , 2]  +  N[3]*DD[I,3]  )  /  (N[I]*DET); 

(*  apply  light-time  correction  and  calculate  time  differences  *) 
RETARD  (JDO, RHO ,  J  D , TAU) ; 

(*  heliocentric  coordinate  vectors  *) 

FOR  I  :=  1  TO  3  DO 
FOR  S  :=  X  TO  Z  DO 

RPL [I . S]  :=  RHO [T] *E [I , S] -RSUN [I , S] ; 

(*  sector/triangle  ratios  eta[i]  *) 

ETA [ 1]  :=  FIND_ETA(  RPL [2] ,  RPL[3],  TAU [ 1 ]  ); 

ETA [2]  :=  FIND_ETA(  RPL [1] .  RPL[3],  TAU[2]  ); 

ETA [3]  :=  FIND_ETA(  RPL [1] ,  RPL [2] .  TAU [ 3 ]  ); 

(*  improvement  of  the  sector/triangle  ratios  *) 

N [1]  :=  (  TAU [1] /ETA [ 1]  )  /  (TAU [2] /ETA [2] ) ; 

N [3]  :=  (  TAU [3] /ETA[3]  )  /  (TAU [2] /ETA [2] ) ; 

WRITELN ( '  rho' , *  1 : 16, RHO [1] : 12: 8, RHO [2] : 12: 8, RHO [3] : 12:8) ; 

UNTIL  (  ABS(RH0[2] -RHOOLD)  <  EPS„RHO  ); 

WRITELN;  WRITELN ( '  Heliocentric  distances  [AU]:');  WRITELN; 

WRITELN  ( ’  r  \*  1  :16, 

NORM (RPL [1]) : 12:8, NORM (RPL [2]) : 12 : 8 , NORM (RPL [3 ] ) : 12 : 8) ; 
WRITELN;  WRITELN; 

(*  derive  orbital  elements  from  first  and  third  observation  *) 
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ELEMENT  ( J D [1] , JD [3] , RPL [1] , RPL [3] ,  TP , Q , ECC , INC , LAN , AOP) ; 


END; 


(* - - - - - - - - - - - *) 

BEGIN 

START (HEADER, RSUN , E , J DO, TEQX) ; 

GAU  SS ( RSUN , E , J  DO , TP , Q , ECC , I NC , LAN , AOP) ; 

DUMPELEM (TP  f Q , ECC , INC , LAN , AOP , TEQX) ; 

SAVEELEM (TP, Q, ECC, INC, LAN, AOP, TEQX, HEADER) ; 

(*  check  solution  *) 

WRITELN ; 

IF  (DOT (E [2] , RSUN [ 2] ) >0)  THEN 

WRITELN  ('  Warning:  observation  in  hemisphere  of  conjunction;', 

1  possible  second  solution'); 

IF  (ECOl.l)  THEN 

WRITELN  ('  Warning:  probably  not  a  realistic  solution  (e>l.l)  '); 

IF  (  (ABS (Q-0.985)<0.1)  AND  (ABS (ECC-0 . 015) <0 . 05)  )  THEN 
WRITELN  ('  Warning:  probably  Earth' 's  orbit  solution'); 

END. 


( 


) 


APPENDIX  G 


Example  ORBDET  Input  File 


1035  Amata:  Dan  Burtz's  Data,  Test  Case  1 
1998  01  21  05.799  03  44  50.43  42  12  41.60  !  Three  observations 
1998  02  13  01.646  03  53  21.41  40  32  49.20  !  Format:  Date  (ymdh), 
1998  03  13  02.213  04  18  34.54  39  20  22.90  !  RA  (hms),  Dec  (hms) 
2000.0  !  Equinox 

2000.0  !  Required  equinox 


APPENDIX  H 


ORBDET  Output  Data 


ORBDET:  orbit  determination  from  three  observations 
version  93/07/01 

(c)  1993  Thomas  Pfleger,  Oliver  Montenbruck 

Summary  of  orbit  determination 
1035  Amata:  Dan  Burtz's  Data,  Test  Case  1 

Initial  data  (ecliptic  geocentric  coordinates  (in  deg)) 

Julian  Date  2450834.74  2450857.57 

Solar  longitude  300.98  324.15 

Planet/Comet  Longitude  63.66  64.91 

Planet/Comet  Latitude  21.81  19.83 

Iteration  of  the  geocentric  distances  rho  [AU] 
rho  2.75015994  3.07983981  3.53209164 

rho  2.68268840  3.01101420  3.44426096 

rho  2.67800252  3.00622284  3.43816414 

rho  2.67766167  3.00587426  3.43772070 

rho  2.67763586  3.00584790  3.43768743 

rho  2.67763398  3.00584598  3.43768499 

rho  2.67763384  3.00584584  3.43768481 

rho  2.67763383  3.00584583  3.43768480 

rho  2.67763383  3.00584582  3.43768480 

Heliocentric  distances  [AU]: 

r  3.28296352  3.32420416  3.37302007 

Orbital  elements  (Equinox  J  2000.00) 

Perihelion  date  tp  1996/ 8/21  13.0623h  (JD  2450317.04) 


Perihelion  distance 

q[AU] 

2.502659 

Semi-major  axis 

a[AU] 

3.139193 

Eccentricity 

e 

0.202770 

Inclination 

i 

18.0870  degrees 

Ascending  node 

Omega 

2.1808  degrees 

Long,  of  perihelion 

Pi 

325.4512  degrees 

Arg.  of  perihelion 

omega 

323.2704  degrees 

Mean  motion 

n 

0.177205  degrees/day 

Mean  anomaly 

M 

85.670976  degrees 

Orbital  period 

P 

5.5621  years 

2450885.59 

352.30 

69.66 

17.67 


ORBDET:  orbit  determination  from  three  observations 
version  93/07/01 

(c)  1993  Thomas  Pfleger,  Oliver  Montenbruck 

Summary  of  orbit  determination 
1035  Amata:  Dan  Burtz's  Data,  Test  Case  2 

Initial  data  (ecliptic  geocentric  coordinates  (in  deg)) 

Julian  Date  2450835.72  2450857.57  2450885.59 

Solar  longitude  301.98  324.15  352.31 

Planet/Comet  Longitude  63.66  64.91  69.66 

Planet/Comet  Latitude  21.73  19.83  17.67 

Iteration  of  the  geocentric  distances  rho  [AU] 

rho  2.76360006  3.08034154  3.53227255 

rho  2.69595418  3.01130373  3.44444420 

rho  2.69126012  3.00650219  3.43835253 

rho  2.69091902  3.00615324  3.43790990 

rho  2.69089249  3.00612615  3.43787602 

rho  2.69089057  3.00612418  3.43787353 

rho  2.69089043  3.00612404  3.43787335 

rho  2.69089042  3.00612403  3.43787333 

rho  2.69089042  3.00612403  3.43787333 

Heliocentric  distances  [AU]: 

r  3.28504090  3.32443219  3.37316824 


Orbital  elements  (Equinox  J  2000.00) 


Perihelion  date 

tp 

1996/8/21  8.8622h  (JD  2450316.87) 

Perihelion  distance 

q[AU] 

2.504413 

Semi-major  axis 

a[AU] 

3.139943 

Eccentricity 

e 

0.202402 

Inclination 

i 

18.0866  degrees 

Ascending  node 

Omega 

2.1676  degrees 

Long,  of  perihelion 

Pi 

325.4894  degrees 

Arg.  of  perihelion 

omega 

323.3218  degrees 

Mean  motion 

n 

0.177142  degrees/day 

Mean  anomaly 

M 

85.671276  degrees 

Orbital  period 

P 

5.5640  years 

ORBDET:  orbit  determination  from  three  observations 
version  93/07/01 

(c)  1993  Thomas  Pfleger,  Oliver  Montenbruck 

Summary  of  orbit  determination 
1035  Amata:  Dan  Burtz's  Data,  Test  Case  3 

Initial  data  (ecliptic  geocentric  coordinates  (in  deg)) 

Julian  Date  2450834.72  2450857.59  2450885.60 

Solar  longitude  300.96  324.17  352.31 

Planet/Comet  Longitude  63.71  64.92  69.66 

Planet/Comet  Latitude  21.29  19.82  17.67 

Iteration  of  the  geocentric  distances  rho  [AU] 

rho  4.19643404  4.05778655  3.93701682 

rho  4.14404303  4.01164478  3.88739874 

rho  4.14216131  4.00998664  3.88561449 

rho  4.14209142  4.00992501  3.88554832 

rho  4.14208882  4.00992272  3.88554586 

rho  4.14208873  4.00992264  3.88554577 

rho  4.14208872  4.00992264  3.88554577 

Heliocentric  distances  [AU]: 
r  4.71512193  4.29416049  3.80401671 


Orbital  elements  (Equinox  J  2000.00) 


Perihelion  date 

tp 

1998/7/18  3.0726h  (JD  2451012.63) 

Perihelion  distance 

q[AU] 

2.489331 

Semi-major  axis 

a[AU] 

-0.733254 

Eccentricity 

e 

4.394909 

Inclination 

i 

18.6437  degrees 

Ascending  node 

Omega 

340.4056  degrees 

Long,  of  perihelion 

Pi 

499.3132  degrees 

Arg.  of  perihelion 

omega 

158.9076  degrees 

Mean  motion 

.  n 

1 .569720  degrees/day 

Mean  anomaly 

M 

-332.981702  degrees 

Orbital  period 

P 

0.6279  years 

Warning:  probably  not  a  realistic  solution  (e>1 .1 ) 


ORBDET:  orbit  determination  from  three  observations 
version  93/07/01 

(c)  1993  Thomas  Pfieger,  Oliver  Montenbruck 

Summary  of  orbit  determination 
1035  Amata:  Dan  Burtz's  Data,  Test  Case  4 

Initial  data  (ecliptic  geocentric  coordinates  (in  deg)) 

Julian  Date  2450840.78  2450857.59  2450885.63 

Solar  longitude  307.13  324.17  352.35 

Planet/Comet  Longitude  63.72  64.92  69.67 

Planet/Comet  Latitude  21.29  19.82  17.67 

Iteration  of  the  geocentric  distances  rho  [AU] 

rho  2.83399704  3.08206213  3.53318200 

rho  2.76532224  3.01193729  3.44535132 

rho  2.76057158  3.00707819  3.43927728 

rho  2.76022749  3.00672620  3.43883734 

rho  2.76020110  3.00669924  3.43880400 

rho  2.76019919  3.00669728  3.43880155 

rho  2.76019905  3.00669714  3.43880138 

rho  2.76019904  3.00669713  3.43880136 

rho  2.76019904  3.00669713  3.43880136 

Heliocentric  distances  [AU]: 

r  3.29445640  3.32471097  3.37356737 


Orbital  elements  (Equinox  J  2000.00) 


Perihelion  date 

tp 

1996/8/21  17.4846h  (JD  2450317.23) 

Perihelion  distance 

q[AU] 

2.504159 

Semi-major  axis 

a[AU] 

3.140728 

Eccentricity 

e 

0.202682 

Inclination 

i 

1 8.0868  degrees 

Ascending  node 

Omega 

2.1664  degrees 

Long,  of  perihelion 

Pi 

325.5467  degrees 

Arg.  of  perihelion 

omega 

323.3803  degrees 

Mean  motion 

n 

0.177076  degrees/day 

Mean  anomaly 

M 

85.575542  degrees 

Orbital  period 

P 

5.5661  years 

ORBDET:  orbit  determination  from  three  observations 
version  93/07/01 

(c)  1993  Thomas  Pfleger,  Oliver  Montenbruck 

Summary  of  orbit  determination 
1035  Amata:  Dan  Burtz's  Data,  Test  Case  5 

Initial  data  (ecliptic  geocentric  coordinates  (in  deg)) 

Julian  Date  2450841.71  2450857.62  2450885.64 

Solar  longitude  308.07  324.20  352.35 

Planet/Comet  Longitude  63.74  64.92  69.67 

Planet/Comet  Latitude  21.20  19.82  17.67 

Iteration  of  the  geocentric  distances  rho  [AU] 

rho  2.84730595  3.08268194  3.53300119 

rho  2.77840944  3.01234947  3.44518341 

rho  2.77364447  3.00747753  3.43911142 

rho  2.77329804  3.00712332  3.43867036 

rho  2.77327278  3.00709749  3.43863820 

rho  2.77327093  3.00709560  3.43863585 

rho  2.77327080  3.00709546  3.43863568 

rho  2.77327079  3.00709545  3.43863567 

rho  2.77327079  3.00709545  3.43863567 

Heliocentric  distances  [AU]: 

r  3.29617835  3.32469276  3.37333148 


Orbital  elements  (Equinox  J  2000.00) 


Perihelion  date 

tP 

1996/8/21  4.97 19h  (JD  2450316.71) 

Perihelion  distance 

q[AU] 

2.506152 

Semi-major  axis 

a[AU] 

3.140653 

Eccentricity 

e 

0.202028 

Inclination 

i 

18.0864  degrees 

Ascending  node 

Omega 

2.1577  degrees 

Long,  of  perihelion 

P> 

325.5288  degrees 

Arg.  of  perihelion 

omega 

323.3711  degrees 

Mean  motion 

n 

0.177082  degrees/day 

Mean  anomaly 

M 

85.670941  degrees 

Orbital  period 

P 

5.5659  years 

ORBDET:  orbit  determination  from  three  observations 
version  93/07/01 

(c)  1993  Thomas  Pfleger,  Oliver  Montenbruck 

Summary  of  orbit  determination 
1035  Amata:  Dan  Burtz's  Data,  Test  Case  6 

Initial  data  (ecliptic  geocentric  coordinates  (in  deg)) 

Julian  Date  2450841.77  2450857.63  2450885.64 

Solar  longitude  308.13  324.21  352.35 

Planet/Comet  Longitude  63.75  64.92  69.67 

Planet/Comet  Latitude  21.20  19.82  17.67 

Iteration  of  the  geocentric  distances  rho  [AU] 

rho  2.84852169  3.08313344  3.53342489 

rho  2.77962613  3.01280417  3.44562781 

rho  2.77486068  3.00793197  3.43955718 

rho  2.77451574  3.00757926  3.43911772 

rho  2.77449272  3.00755568  3.43908782 

rho  2.77449102  3.00755394  3.43908565 

rho  2.77449090  3.00755382  3.43908550 

rho  2.77449089  3.00755381  3.43908549 

Heliocentric  distances  [AU]: 

r  3.29668444  3.32509755  3.37373549 


Orbital  elements  (Equinox  J  2000.00) 


Perihelion  date 

<P 

1996/8/21  7.91 03h  (JD  2450316.83) 

Perihelion  distance 

q[AU] 

2.507582 

Semi-major  axis 

a[AU] 

3.141972 

Eccentricity 

e 

0.201908 

Inclination 

i 

18.0861  degrees 

Ascending  node 

Omega 

2.1416  degrees 

Long,  of  perihelion 

P' 

325.6097  degrees 

Arg.  of  perihelion 

omega 

323.4682  degrees 

Mean  motion 

n 

0.176970  degrees/day 

Mean  anomaly 

M 

85.595332  degrees 

Orbital  period 

P 

5.5694  years 

APPENDIX  I 


FIND  ORB  Input  File 


J98X01X  Cl  998  01  21.24164  03  44 
J98X01X  Cl  998  01  21.24550  03  44 
J98X01X  Cl  998  01  21 .27009  03  44 
J98X01X  C1998  01  21.27427  03  44 
J98X01X  Cl  998  01  22.21771  03  44 
J98X01X  Cl  998  01  22.22047  03  44 
J98X01X  C1998  01  22.27830  03  44 
J98X01X  C1998  01  22.28038  03  44 
J98X01X  C1998  01  27.21590  03  45 
J98X01X  Cl  998  01  27.21848  03  45 
J 98X0 IX  C1998  01  27.25006  03  45 
J98X01X  Cl  998  01  27.25433  03  45 
J98X01X  C1998  01  27.27655  03  45 
J98X01X  C1998  01  27.28177  03  45 
J98X01X  Cl  998  01  28.21384  03  46 
J98X01X  C1998  01  28.21614  03  46 
J98X01X  Cl  998  01  28.24293  03  46 
J98X01X  Cl  998  01  28.24633  03  46 
J98X01X  Cl  998  01  28.26623  03  46 
J98X01X  C1998  01  28.27069  03  46 
J98X01X  Cl  998  02  13.06861  03  53 
J98X01X  C1998  02  13.07160  03  53 
J98X01X  Cl  998  02  13.09028  03  53 
J98X01X  Cl 998  02  13.09264  03  53 
J98X01X  Cl  998  02  13.12376  03  53 
J98X01X  Cl  998  02  13.12638  03  53 
J98X01X  Cl  998  03  13.09222  04  18 
J98X01X  Cl  998  03  13.09472  04  18 
J 98X0 IX  C1998  03  13.09684  04  18 
J98X01X  C1998  03  13.13267  04  18 
J98X01X  C1998  03  13.13863  04  18 
J98X01X  Cl  998  03  13.14065  04  18 


50.43 

+42 

12 

41.6 

15.9 

V 

712 

50.41 

+42 

12 

40.7 

15.9 

V 

712 

50.55 

+42 

12 

32.5 

15.9 

V 

712 

50.54 

+42 

12 

31.2 

16.0 

V 

712 

55.41 

+42 

07 

42.6 

15.9 

V 

712 

55.42 

+42 

07 

41.2 

15.9 

V 

712 

55.73 

+42 

07 

23.9 

16.1 

V 

712 

55.74 

+42 

07 

23.1 

16.1 

V 

712 

46.10 

+41 

43 

02.5 

16.5 

V 

712 

46.10 

+41 

43 

01.6 

16.6 

V 

712 

46.54 

+41 

42 

53.6 

15.7 

V 

712 

46.62 

+41 

42 

51.9 

16.0 

V 

712 

46.94 

+41 

42 

45.1 

16.0 

V 

712 

47.03 

+41 

42 

43.5 

16.0 

V 

712 

01.17 

+41 

38 

18.7 

16.1 

V 

712 

01.20 

+41 

38 

18.3 

15.9 

V 

712 

01.63 

+41 

38 

10.2 

16.0 

V 

712 

01.68 

+41 

38 

09.1 

16.1 

V 

712 

01.96 

+41 

38 

03.8 

16.0 

V 

712 

02.02 

+41 

38 

02.5 

16.0 

V 

712 

21.41 

+40 

32 

49.2 

16.6 

V 

712 

21.53 

+40 

32 

48.7 

16.5 

V 

712 

22.22 

+40 

32 

44.8 

16.4 

V 

712 

22.33 

+40 

32 

44.0 

16.4 

V 

712 

23.49 

+40 

32 

37.6 

16.4 

V 

712 

23.59 

+40 

32 

37.1 

16.4 

V 

712 

34.54 

+39 

20 

22.9 

16.2 

V 

712 

34.71 

+39 

20 

22.4 

16.4 

V 

712 

34.84 

+39 

20 

22.2 

16.1 

V 

712 

37.23 

+39 

20 

18.7 

16.2 

V 

712 

37.63 

+39 

20 

17.6 

16.3 

V 

712 

37.76 

+39 

20 

17.3 

16.1 

V 

712 

APPENDIX  J 


FIND  ORB  Output  File 


Orbital  elements: 

1998XX1 

Perihelion  1996  Aug  21.102725  TT 
Epoch  1997  Dec  18.0  TT  =  JDT  2450800.5 
M  85.82430  (2000.0)  P  Q 

n  0.17736059  Peri.  323.13457  0.82135348  0.57029497 

a  3.1373621  Node  2.20044  -0.42057277  0.61957291 

e  0.2025591  Incl.  18.08732  -0.38535309  0.53934502 

P  5.56  H  10.4  G  0.15  q  2.5018609 
From  32  observations  1998  Jan.  21 -Mar.  13;  RMS  error  0.360  arcseconds 


APPENDIX  K 


Example  Astrometry  Submission  to  MPC 


USAF  Academy  Observatory 

1998  January  26 

Observatory  Code:  712 

Longitude:  1O4.88110W  p  cos  i': +0.77837 

Latitude:  39.OO7O0N  p  sin  i':  +0.62625 

Height:  2180.0m 

41 -cm  reflector  +  CCD,  GSC  reference  stars  (J2000.0). 
Postal  adress: 

Wetterer 

Department  of  Physics 
USAF  Academy,  Colorado 
USA 


Object 

Date 

U.T. 

R.A. 

(2000.0)  Decl. 

Mag. 

Observer  Comp  Stars 

Karayusuf 

1998 

01 

24.41848 

10 

13 

41.85 

+35 

13 

28.8 

16.1V 

Wetterer 

8 

Karayusuf 

1998 

01 

24.42039 

10 

13 

41.75 

+35 

13 

32.6 

15.9V 

Wetterer 

8 

Karayusuf 

1998 

01 

24.43358 

10 

13 

41.27 

+35 

13 

53.9 

15.9V 

Wetterer 

8 

Karayusuf 

1998 

01 

24.44428 

10 

13 

40.99 

+35 

14 

09.5 

15.8V 

Wetterer 

8 

Amata 

1998 

01 

27.21590 

03 

45 

46.10 

+41 

43 

02.5 

16.5V 

Wetterer/Burtz 

9 

Amata 

1998 

01 

27.21848 

03 

45 

46.10 

+41 

43 

01.6 

16.6V 

Wetterer/Burtz 

9 

Amata 

1998 

01 

27.25006 

03 

45 

46.54 

+41 

42 

53.6 

15.7V 

Wetterer/Burtz 

9 

Amata 

1998 

01 

27.25433 

03 

45 

46.62 

+41 

42 

51.9 

16.0V 

Wetterer/Burtz 

9 

Amata 

1998 

01 

27.27655 

03 

45 

46.94 

+41 

42 

45.1 

16.0V 

Wetterer/Burtz 

8 

Amata 

1998 

01 

27.28177 

03 

45 

47.03 

+41 

42 

43.5 

16.0V 

Wetterer/Burtz 

8 

Amata 

1998 

01 

28.21384 

03 

46 

01.17 

+41 

38 

18.7 

16.1V 

Burtz/Wetterer 

9 

Amata 

1998 

01 

28.21614 

03 

46 

01.20 

+41 

38 

18.3 

15.9V 

Burtz/Wetterer 

9 

Amata 

1998 

01 

28.24293 

03 

46 

01.63 

+41 

38 

10.2 

16.0V 

Burtz/Wetterer 

9 

Amata 

1998 

01 

28.24633 

03 

46 

01.68 

+41 

38 

09.1 

16.1V 

Burtz/Wetterer 

9 

Amata 

1998 

01 

28.26623 

03 

46 

01.96 

+41 

38 

03.8 

16.0V 

Burtz/Wetterer 

9 

Amata 

1998 

01 

28.27069 

03 

46 

02.02 

+41 

38 

02.5 

16.0  V 

Burtz/Wetterer 

9 

Ramanujan 

1998 

01 

27.22343 

05 

26 

21.20 

+10 

25 

24.6 

17.3V 

Wetterer/Burtz 

12 

Ramanujan 

1998 

01 

27.28451 

05 

26 

20.13 

+10 

25 

40.3 

17.7V 

Wetterer/Burtz 

12 

Ramanujan 

1998 

01 

27.28787 

05 

26 

19.97 

+10 

25 

39.8 

17.5V 

Wetterer/Burtz 

12 

Ramanujan 

1998 

01 

28.22064 

05 

26 

03.15 

+10 

28 

49.1 

17.4V 

Burtz/Wetterer 

12 

Ramanujan 

1998 

01 

28.22256 

05 

26 

03.08 

+10 

28 

49.1 

17.4V 

Burtz/Wetterer 

12 

Ramanujan 

1998 

01 

28.24888 

05 

26 

02.53 

+10 

28 

54.3 

17.1V 

Burtz/Wetterer 

12 

Ramanujan 

1998 

01 

28.25224 

05 

26 

02.53 

+10 

28 

55.0 

17.3V 

Burtz/Wetterer 

12 

Comrie 

1998 

01 

27.22972 

10 

04 

52.35 

+19 

04 

12.7 

17.4V 

Wetterer/Burtz 

10 

Comrie 

1998 

01 

27.23221 

10 

04 

52.22 

+19 

04 

14.1 

17.1V 

Wetterer/Burtz 

10 

Comrie 

1998 

01 

27.26237 

10 

04 

50.62 

+19 

04 

25.4 

17.0V 

Wetterer/Burtz 

10 

Comrie 

1998 

01 

27.26634 

10 

04 

50.43 

+19 

04 

27.4 

17.0V 

Wetterer/Burtz 

10 

Comrie 

1998 

01 

28.22870 

10 

04 

01.50 

+19 

10 

32.9 

17.1V 

Burtz/Wetterer 

8 

Comrie 

1998 

01 

28.23061 

10 

04 

01.32 

+19 

10 

34.2 

17.2V 

Burtz/Wetterer 

8 

Comrie 

1998 

01 

28.25449 

10 

04 

00.03 

+19 

10 

43.8 

16.9V 

Burtz/Wetterer 

8 

Comrie 

1998 

01 

28.25843 

10 

03 

59.82 

+19 

10 

45.3 

17.0V 

Burtz/Wetterer 

8 

Karayusaf  1998  01  27.23680  10  12  01 

Karayusaf  1998  01  27.23909  10  12  01 

Karayusaf  1 998  01  27.26940  1 0  1 2  00 

Karayusaf  1998  01  27.27315  10  12  00 

Karayusaf  1998  01  28.23657  10  11  19 

Karayusaf  1998  01  28.23846  10  11  19 

Karauysaf  1998  01  28.26065  10  11  18 

Karauysaf  1998  01  28.26385  10  11  18 

Raffinetti  1998  01  27.29867  11  09  48 

Raffinetti  1998  01  27.30102  11  09  48 

Raffinetti  1998  01  27.33069  1 1  09  48 

Raffinetti  1998  01  27.33454  11  09  48 

Raffinetti  1998  01  28.29027  11  09  32 

Raffinetti  1998  01  28.30932  11  09  32 

Raffinetti  1998  01  28.31381  11  09  31 

Pien  1998  01  27.30875  11  22  51 

Pien  1998  01  27.31174  11  22  51 

Pien  1998  01  27.33708  1 1  22  50 

Pien  1998  01  27.34104  11  22  50 

Pien  1998  01  28.27936  11  22  28 

Pien  1998  01  28.28119  11  22  28 

Pien  1998  01  28.30284  11  22  27 

Pien  1998  01  28.30668  11  22  27 

Pien  1998  01  28.33435  1 1  22  27 

Pien  1998  01  28.32125  1 1  22  27 


.87+36  32  02.6  15.4V  Wetterer/Burtz  9 
.75+36  32  06.6  15.4V  Wetterer/Burtz  9 
.38+36  32  57.8  15.4V  Wetterer/Burtz  9 
.28+36  33  03.1  15.4V  Wetterer/Burtz  9 
.86+36  59  40.8  16.0V  Burtz/Wetterer  7 
.76+36  59  43.7  16.0V  Burtz/Wetterer  7 
.71  +37  00  21.0  15.9V  Burtz/Wetterer  7 
.55+37  00  26.3  16.0V  Burtz/Wetterer  7 
.68-05  20  32.3  16.6V  Wetterer/Burtz  10 
.59-05  20  32.0  16.6V  Wetterer/Burtz  10 
.09-05  20  33.2  16.8V  Wetterer/Burtz  10 
.08-05  20  32.7  16.9V  Wetterer/Burtz  10 
.38-05  20  26.6  16.8V  Burtz/Wetterer  10 
.00-05  20  26.6  16.8V  Burtz/Wetterer  10 
.92-05  20  26.3  16.9V  Burtz/Wetterer  10 
.58+13  54  15.8  17.5V  Wetterer/Burtz  6 
.58  +13  54  15.9  17.2V  Wetterer/Burtz  6 
.94+13  54  26.7  17.5V  Wetterer/Burtz  6 
.67+13  54  26.1  16.9V  Wetterer/Burtz  6 
.55+13  59  47.9  16.8V  Burtz/Wetterer  7 
.51  +13  59  48.1  16.9V  Burtz/Wetterer  7 
.92+13  59  56.3  16.5V  Burtz/Wetterer  7 
.87+13  59  57.5  16.6V  Burtz/Wetterer  7 
.15+14  00  07.0  16.5V  Burtz/Wetterer  7 
.51  +14  00  02.2  16.5V  Burtz/Wetterer  7 


